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Many biological systems are modulated by unknown slow processes. This can severely 
hinder analysis - especially in excitable neurons, which are highly non-linear and stochastic 
systems. We show the analysis simplifies considerably if the input matches the sparse 
"spiky" nature of the output. In this case, a linearized spiking Input-Output (I/O) relation 
can be derived semi-analytically, relating input spike trains to output spikes based on 
known biophysical properties. Using this I/O relation we obtain closed-form expressions 
for all second order statistics (input - internal state - output correlations and spectra), 
construct optimal linear estimators for the neuronal response and internal state and 
perform parameter identification. These results are guaranteed to hold, for a general 
stochastic biophysical neuron model, with only a few assumptions (mainly, timescale 
separation). We numerically test the resulting expressions for various models, and show 
that they hold well, even in cases where our assumptions fail to hold. In a companion 
paper we demonstrate how this approach enables us to fit a biophysical neuron model so 
it reproduces experimentally observed temporal firing statistics on days-long experiments. 



Keywords: conductance based neuron models, noise, ion channels, adaptation, power spectral density, linear 
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1. INTRODUCTION 

Neurons are modeled biophysically using Conductance-Based 
Models (CBMs). In CBMs, the membrane time constant and 
the timescales of fast channel kinetics determine the timescale 
of Action Potential (AP) generation in the neuron. These are 
typically around 1-20 ms. However, there are various modulat- 
ing processes that affect the response on slower timescales. Many 
types of ion channels exist, and some change with a timescale 
as slow as 10 s (Channelpedia), and possibly even minutes (Toib 
et al., 1998). Additional new sub-cellular kinetic processes are 
being discovered at an explosive rate (Bean, 2007; Sjostrom et al., 
2008; Debanne et al., 2011). This variety is particularly large for 
very slow processes (Marom, 2010). 

For example, ion channels are known to be regulated over 
the course of long timescales (Levitan, 1994; Staub et al., 1997; 
Jugloff, 2000; Monjaraz et al., 2000), which could cause changes 
in ion channel numbers, conductances and kinetics. Also, the 
ionic concentrations in the cell depend on the activity of the ionic 
pumps, which can be affected by the metabolism of the network 
(Silver et al, 1997; Kasischke et al, 2004). Finally, the cellular neu- 
rites (De Paola et al., 2006; Nishiyama et al., 2007) and even the 
spike initiation region (Grubb and Burrone, 2010) can shift their 
location with time. All these changes can have a large effect on 
excitability. 

Therefore, current CBMs can be considered as strictly accurate 
only below a certain timescale, since they do not incorporate most 
of these slow processes. A main reason for this "neglect" is that 
such slow processes are not well characterized. This is especially 
problematic since neurons are excitable, so their dynamics is far 
from equilibrium, highly non-linear and contain feedback. Due 



to the large number of processes which are unknown or lacking 
known parameters, it would be hard to simulate or analyze such 
models. Therefore, it may be hard to quantitatively predict how 
adding and tuning slow processes in the model would affect the 
dynamics at longer timescales. 

In order to allow CBMs with many slow process to be fit- 
ted and analyzed, it is desirable to have general expressions that 
describe their Input-Output (I/O) relation explicitly, based on 
biophysical parameters. In a previous paper (Soudry and Meir, 
2012b), we found that this becomes possible if we use (experi- 
mentally relevant Elul and Adey, 1966; Kaplan et al., 1996; De Col 
et al., 2008; Gal et al, 2010; Goldwyn et al, 2012) sparse spike 
inputs, similar to the typical output of the neuron (Figures 1A,B). 
In this case, we derived semi-analytically 1 a discrete piecewise 
linear map describing the neuronal dynamics between stimula- 
tion spikes, for a general deterministic neuron model with a few 
assumptions (mainly, a timescale separation assumption). Based 
on this reduced map, we were able to derive expressions for the 
"mean" behavior of the neuron (e.g., firing modes, firing rate and 
mean latency). 

In this paper, we find that stronger and more general analytical 
results can be obtained if we take into account the stochastic - 
ity of the neuron - arising from ion channel noise 2 (Neher 
and Sakmann, 1976; Hille, 2001). Due to the presence of this 
noise, the discrete map describing the neuronal dynamics is 



A semi-analytic derivation is an analytic derivation in which some terms are 
obtained by relatively simple numerics. See 2.2 for our implementation. 
2 We demonstrated that such noise should strongly affect the neuronal 
response to sparse stimulation (Soudry and Meir, 2012b). 
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FIGURE 1 | Schematic summary. (A) Aim: find the I/O relation 
between inter-stimulus intervals (T m ) and Action Potential (AP) 
occurrences (V m ) - for a general biophysical neuron model 
(Equations 1-3). (B) An AP "occurred" if the voltage V crossed 



threshold V th following the (sparse) stimulus, with 7 m 2>r AP . (C) 
Result: Biophysical neuron model reduced to a simple linear system 
with feedback (Equations 11, 12), and biophysically meaningful 
parameters (F,d,a and w). 



"smoothed out," and can be linearized. This linearized map con- 
stitutes a concise description for the neuronal I/O (Equations 
11, 12) based on biophysically meaningful parameters. This I/O 
is well described by an "engineering-style" block diagram with 
feedback (Figure 1C), where the input is the process of stimu- 
lation intervals and the output is the AP response (Figure 1A). 
Note that the response is affected both by internal noise and by 
the input. Beyond conceptual lucidity, such a linear I/O allows 
the utilization of well known statistical tools to derive all sec- 
ond order statistics, to construct linear optimal estimators and to 
perform parameter identification. These results hold numerically 
(Figure 3), even sometimes when our assumptions break down 
(Figure 4). 

In our previous paper, Soudry and Meir (2012b), we used 
our results to model recent experiments (Gal et al., 2010) where 
synaptically isolated individual neurons, from rat cortical cul- 
ture, were stimulated with extra-cellular sparse current pulses for 
an unprecedented duration of days. Our results enabled us to 
explain the "mean" response of these neurons. However, the sec- 
ond order-statistics in the experiment seem particularly puzzling. 
The neurons exhibited 1//" statistics (Keshner, 1982), respond- 
ing in a complex and irregular manner from seconds to days. In a 
companion paper (Soudry and Meir, 2014), we demonstrate the 
utility of our new results. These results allow us to reproduce and 
analyze the origins of this 1 /f* on very long timescales. 

2. RESULTS 

This section described our main results in outline. The details of 
each sub-section here appear in the corresponding sub-section 
in section 4. For readers who do not wish to go through the 
detailed derivations, the present section is self-contained. Readers 
who do wish to follow the mathematical derivations, should first 



read section 4, where, for convenience, each subsection (except 
for the last one) can be read independently. In our notation (•> 

is an ensemble average, i = *J — 1, a non-capital boldfaced let- 
ter x = (xi, ... , x n ) T is a column vector (where (-) T denotes 
transpose), and a boldfaced capital letter X is a matrix (with 
components X mn ) . 

2.1. FULL MODEL 

The voltage dynamics of an isopotential neuron are determined 
by ion channels, protein pores which change their conforma- 
tions stochastically with voltage-dependent rates (Hille, 2001). At 
the population level, such dynamics are generically described by 
Fox and Lu (1994), Goldwyn et al. (2011), and Soudry and Meir 
(2012b) a CBM 

V=f(V,t, S ,I(t)) (1) 
r = A r (V)r + B r (V,r)? r (2) 
s = A s (V)s + B s <y,s)§ s , (3) 

with voltage V, stimulation current I(t), rapid variables r 
(e.g., m, n, h in the Hodgkin-Huxley (HH) model Hodgkin and 
Huxley, 1952), slow "excitability" variables s (e.g., slow sodium 
inactivation Chandler and Meves, 1970), rate matrices A r / S , white 
noise processes £ r / s (with zero mean and unit variance), and 
matrices B r / S which can be written explicitly using the rates 
and ion channel numbers (Orio and Soudry, 2012) (D = BB T 
is the diffusion matrix Orio and Soudry, 2012). For simplicity, 
we assumed that r and s are not coupled directly, but this is 
non-essential (Contou-Carrere, 2011; Wainrib et al, 2011). The 
parameter space can be constrained (Soudry and Meir, 2012b), 
since we consider here only excitable, non-oscillatory neurons 
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which do not fire spontaneously 3 and which have a single resting 
state - as is common for cortical cells, e.g., Gal et al. (2010). 

Since the components of r and s usually represent fractions, 
in some cases it is more convenient to use the normalization con- 
straint (i.e., that fractions sum to one), and reduce the dimensions 
of r, s, and f r / s . After this reduction, the form of Equations (1-3) 
changes to 

V =f(V,r,s,I(t)) (4) 
r = A r (V)r-b r (V)+B r (V,r)§ r , (5) 
s = A s (V) s - b 5 {V) + B 5 (V, s) | s , (6) 

where all the variables and parameters have been redefined (with 
their size decreased). Note that we have slightly abused nota- 
tion by using the same symbols in Equations (4-6) and in 
Equations (1-3). The specific set of equations used will always 
be stated. We call Equations (4-6) the "compressed form" of 
the CBM. 

Such biophysical neuronal models (either Equations 1-3 or 
4-6) are generally complex and non-linear, containing many 
variables and unknown parameters (sometimes ranging in the 
hundreds Koch and Segev, 1989; Roth and Hausser, 2001), not 
all of which can be identified (Huys et al., 2006). Therefore, such 
models are notoriously difficult to tune, highly susceptible to 
over-fitting and computationally expensive (Migliore et al., 2006; 
Gerstner and Naud, 2009; Druckmann et al., 201 1). Also, the high 
degree of non-linearity usually prevents exact mathematical anal- 
ysis of such models at their full level of complexity (Ermen trout 
and Terman, 2010). However, much of the complexity in such 
models can be overcome under well defined and experimentally 
relevant settings (Elul and Adey, 1966; Kaplan et al, 1996; De Col 
et al., 2008; Gal et al, 2010; Goldwyn et al, 2012), if we use sparse 
inputs, similar in nature to the spikes commonly produced by the 
neuron. 

2.2. MODEL REDUCTION 

We consider a stimulation setting motivated by the experiments 
described in Gal et al. (2010) and further elaborated on in sec- 
tion 3. Specifically, suppose I (r) consists of a train of pulses 
arriving at times {t m } (Figure 1A, top), so T m = t m+ \ — t m ^> 
tap with tap being the timescale of an AP (Figure IB). Our aim 
is to describe the AP occurrences Y m , where Y m = 1 if an AP 
occurred immediately after the m-th stimulation, and 0 other- 
wise (Figure 1A, bottom). Recall again that we assume the neuron 
does not generate APs unless stimulated (as observed in Gal et al., 
2010). 

In this section we "average out" Equations (1-3) using a semi- 
analytical method similar to that in Soudry and Meir (2012b).To 
do so, we need to integrate Equations (1-3) between r m and 
t m +i. Since T m ^> tap, the rapid AP generation dynamics of 
(V, r) relax to a steady state before t m + \. Therefore, the neuron 
AP "remembers" any history before t m only through s m = s (t m ). 
Given s m , the response of the fast variables (V, r) to the m-th 



3 I.e., if Vt : I (t) = 0, then the probability that a neuron will fire is negligible - 
on any relevant finite time interval (e.g., minutes or days). 



stimulation spike will determine the probability to generate an 
AP. This probability, 

Pap (Sm) = P (Y m \s m ) = (Y m \s m ) , 

collapses all the relevant information from Equations (1,2), and 
can be found numerically from the pulse response of Equations 
(1, 2) with s held fixed (section 4.2.4). 

In order to integrate the remaining Equation (3), we define 
X + , X_ and Xq to be the averages of a quantity X 5 during an AP 
response, a failed AP response and rest, respectively 4 . Also, we 
denote 

X(Y m , T m ) = tap^ 1 {Y m X+ + (1 - F m )X_) + 

(1 - TapT- 1 )^, (7) 

as the steady state mean value of X s . For analytical simplicity we 
assume 5 T m <SC t s . We obtain, to first order 

®m+ 1 = Sm ~t~ T m A (Y m , T m ) s m -\- n m . (8) 

where n m is a white noise process with zero mean and vari- 
ance T m D(Y m , T m ). For the compressed form (Equations 4-6) 
we have instead 

Sm+ l = &m "f" T m ^A (Y m , T m ) s m b (Y m , T m )J -j- n m . (9) 

Note that such a simplified discrete time map, which describes the 
excitability dynamics of the neuron, has far fewer parameters than 
the full model, since it is written explicitly only using the aver- 
aged microscopic rates of s (through A and D), population sizes 
(through D), the probability to generate an AP given s, pAP (s), 
and the relevant timescales. This effective model exposes the large 
degeneracy in the parameters of the full model and leads to signif- 
icantly reduced simulation times and mathematical tractability. 
Notably, the dynamics of the state s m (Equation 8) depends on 
the input T m and the output Y m - and this feedback affects all of 
our following results. 

2.3. LINEARIZATION 

In this section we exploit the intrinsic ion channel noise to 
linearize the neuronal dynamics, rendering it more tractable 
than the (less realistic) noiseless case (Soudry and Meir, 2012b). 
Suppose that the inter-stimulus intervals {T m } have stationary 
statistics with mean T* so that tap <C T m t s with high prob- 
ability. Since s is slow and AP generation is rather noisy in this 
regime (Soudry and Meir, 2012b) (so pAP (s m ) is slowly varying), 
we assume that a stable excitability fixed point s* exists (Figure 2). 
Therefore, the perturbations s m = s m — s* are small and we can 
linearize 

PAP (Sm) ~ P* + W T S m . 



4 E.g., as in Equations (50-52). Note also a similar notation was also used in 
Soudry and Meir (2012b) (e.g., Equations 2.15, 2.16), where we used H/M/L 
instead of +/ — /0. 

5 Later we shall demonstrate numerically that this is not a necessary condition. 
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FIGURE 2 | Schematic explanation of linearization. In a deterministic 
neuron, an AP will be generated in response to stimulation if and only if the 
neuronal excitability (here, s) is above a certain threshold (A). This 
generates discontinuous dynamics in the neuronal excitability (B), see 
Equations 7, 8). In a stochastic neuron, the response probability is a smooth 
function of s (C). In turn, this "smooths" the dynamics (D). Note that if the 
noise is sufficiently high (as is true in many cases, for biophysically realistic 
levels of noise), then this generates a stable fixed point s* - which gives 
the mean response probability p„, and around which the dynamics can be 
linearized (yellow region). 



Denoting X* = X (p*, T*), the mean AP firing rate can be found 
self consistently from the location of the fixed point s*, 

{Y m ) = p* = Pap (s*) , (10) 

where s* depends on p* through A^s* = 0, or s* = A^'b* in the 
compressed form. 

The perturbations s m = s m — s* around the fixed point s* are 
described by the linear system 

Sm + 1 = Fs m + dT m + aY m + n M , (11) 
Y m = w T s m + e m , (12) 

where f,„ = T m - X*, Y m = Y m - (Y m ), F = I + T^A*, 
(n m n„,) = T^D*, e m is a (non-Gaussian) white noise process, 

(e m ) = (e m n m ) = 0, <r 2 = (e 2 m ) = p* (l - p*), d = Ans* and 

a = tap (A+ — A_) s*. If we use the compressed form instead, 
then these results remain valid, except we need to re-define 

d = Ans* - b 0 and a = t A p [(A+ - A_) s* - (b + - b_)]. 

The linear I/O for the fluctuations in Equations (11, 12), 
which contains feedback from the "output" Y m to the state vari- 
able s m (Figure 1C), can be very helpful mathematically and its 
parameters are directly related to biophysical quantities. 

2.4. LINEAR SYSTEMS ANALYSIS 

Using standard tools, this formulation makes it now possible to 
construct optimal linear estimators for Y m and s m (Anderson 
and Moore, 1979), perform parameter identification (Lejung, 
1999), and find all second order statistics in the system (Papoulis 
and Pillai, 1965; Gardiner, 2004), such as correlations or Power 



Spectral Densities (PSD). For example, for / 7^ , the PSD of 
the output is 

S Y (f) = w T H c (-/) (D* + T- 2 dd T S r (f)) Hj (f) w (13) 

+ r„ ( r c 2 |i + r- 1 w T H c (f)a| 2 

where 

H, (f) ± (litfi - A, - T-W 7 )" 1 . (14) 

Similarly, the PSD of the state variables is 

S s {f) = H c (-/) (D, + T-'aaTa 2 + T- 2 dd T S r (fj) Hj (f) , 

(15) 

and the input-output cross-PSD is 

Syr if) = X- x w T H c (-/) dS T (/) . (16) 

Again, note the large degeneracy here - many different sets of 
parameters will generate the same PSD. Using similar methods, 
the PSDs of various response features, such as the AP latency or 
amplitude, can also be derived (Equation 124). 

Finally, we note Equations (11) and (12) can be re-arranged as 
a direct I/O relation. First, we define the filters (transfer functions) 

H ext (f) 4 I^VH^d (17) 
H int if) = (x-'w T H c (/) K + l) a v (18) 

where K = a + FPwcr^ 2 and er 2 = w T Pw + c 2 ,, with P being 
the solution of 

P = FPF T - (w T Pw + cr 2 ) 1 FPww T PF T + T^D*. (19) 

Using these filters, we obtain, in the frequency domain, 

Y{f)=H^(f)T(f)+H int (f)z(f), (20) 

where Y (f ), T (f) and z(f) are the Fourier transforms of Y m , T m 
and z m , respectively, with z m being a white noise process with 
zero mean and unit variance. Notably, these transfer functions 
can be identified from the spiking input-output of the neuron 

| T m , Y m J , without access to the underlying dynamics or biophys- 
ical parameters. Specifically, Equation (20) has the form of an 
ARMAx(M, M, M) model 6 (Lejung, 1999) (recall M is the dimen- 
sion of s), which can be estimated using standard tools (e.g., the 
system identification toolbox in Matlab). 



6 Note a more general Box-Jenkins model is not required, since the poles of 
H ext (/) and H mt (/) are identical (assuming no pole-zero cancelation). 
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2.5. NUMERICAL TESTS 

As we argued so far, a main asset of the present approach is its 
applicability to a broad range of models of various degrees of 
complexity and realism. Recall that our three assumptions are 

1. tap <3C T m (temporally sparse input). 

2. T m <SC f s (timescale separation). 

3. A stable excitability fixed point s* exists, ("noisy" neuron). 

In this section we will demonstrate that our analytical approxi- 
mations agree very well with the numerical solution of Equations 
(1-3), even in some cases where the assumptions 2 and 3 do 
not hold. Therefore, these assumptions are sufficient, but not 
necessary. 

2.5.1. The HHS model 

First, in Figure 3 we tested our results on the HH model with 
Slow sodium inactivation. This "HHS" model (Soudry and Meir, 
2012b, and see section 4.5.1 for parameter values) augments the 
classic HH model (Hodgkin and Huxley, 1952) with an additional 
slow inactivation process of the sodium conductance (Chandler 
and Meves, 1970; Fleidervish et al., 1996). The HHS model 
includes the uncoupled stochastic Hodgkin-Huxley (HH) model 
equations (Fox and Lu, 1994), and is written in the compressed 



formulation (Equations 4-6) 

CV = g Na sm 3 h (E Na -V)+ g K n 4 (E K - V) 

+ gL(E L -V)+I(t) (21) 
r = [ctr {V) (1 - r) - 0, (V) r]4> + 

VN-V (a r (V) (1 - r) + J r (V) r)f r , (22) 

for r = m, n and h, with the additional kinetic equation for slow 
sodium inactivation 

s = S (V) (1 -s) - y (V)s + JlST 1 {8 (V) (1 - s) + y (V)s)&, 

(23) 

where V is the membrane voltage, I (f) is the input current, m, n 
and h are ion channel "gating variables," a r (V) , f) r (V) , S (V) , 
and y (V) are the voltage dependent kinetic rates of these gating 
variables, (p is an auxiliary dimensionless number, C is the mem- 
brane's capacitance, Ek, E^ a and Ei are ionic reversal potentials, 
gK, gNa and gi are ionic conductances and N is the number of ion 
channels. Note that in this model r s is between 20 s (at rest) and 
40 s (during an AP). 

In Figure 3A we show that through Equation (10) we can 
accurately calculate , the mean probability to generate an AP 
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FIGURE 3 | Comparing the mathematical results with the numerical 
simulation of the full model (Equations 1-3) for the stochastic HHS 
model (section 4.5.1). (A) Firing probability p» (T" 1 ) (Equation 10) for 
different currents (/ stim = 7.5, 7.7, 7.9, 8.1. 8.3 /jt,A from bottom to top). (B) 
The PSDs S Y (f) and S s (f). "Sim" is a simulation of the full model, "Map" is 
a (10 4 faster) simulation of Equation (8) together with pap (s m ), while 
"Approx" refers to the analytical expressions (Equations 13-15). "Ident" is 



the PSD Sy (f) of the linear system identified from the spiking data. Note the 
high/low-pass filter shapes of Sy(f) and S s (f), respectively. (C) Optimal linear 
estimation of s. (D) Amplitude and phase of the cross-spectrum Syr {f) for 
Poisson stimulation (Equations 16). Note that the frequency range was cut 
due to spectral estimation noise (see Figure 8). Parameters: Iq = 7.9 ij,A and 
T» = 50 ms in (B-D), and also stimulation is periodical in (A-C). Note the 
low-pass filter shapes of SYr( f)- 
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(so p*!^ 1 is the firing rate of the neuron). In Figure 3B we 
demonstrate both the analytical expression (Equations 13, 15), 
or a simulation of the reduced model (Equation 8), will give 
the PSDs Sy (f) or S 5 (/) of the full model (Equations 1-3). In 
Figure 3D we do the same for the analytical expression (Equation 
16) of the Cross-PSD Syr (/)• In Figure 3C we show that we can 
construct a linear optimal filter for the internal state s m , given 

I { Tjtj^Tp 1 , {yjtl^To 1 J quite well, with low mean square error 
(section 4.4.4). Finally, back in Figure 3B, top, we infer the linear 
model parameters from the spike output using system identifi- 
cation tools [here, with ARMAx(l, 1, 1)], and present the PSD 
of the identified model ("Ident"). Since Sy (f) = \Hi nt (/) 2 (see 
Equation 111) for periodical input (in which T m = 0) this allows 
us to confirm that the linear model was identified. As can be seen, 
the identified filter matches well with that of the linear system. 

2.5.2. Testing the limit of our assumptions 

Next, we demonstrate that our analytical expressions hold also for 
various other models. Specifically, in the following scenarios: (1) 
when the kinetics of the neuron are extended to arbitrarily slow 
timescales, (2) when the assumptions 2 and 3 break down, (3) 
when the rapid and slow kinetics are coupled, (4) when "physi- 
ological" synaptic inputs are used. These results are presented in 
Figures 4, 5, with specific model parameters given in section 4.5. 
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First, we tested whether or not the model can be extended 
to arbitrarily slow timescales. We added to the HHS model four 
types of slow sodium inactivation processes with increasingly 
slower kinetics and smaller channel numbers. In the first case, 
those processes were added additively (as different currents), so 
s was replaced with ^ s; in the voltage equation (Equation 21). 
This model was denoted "HHMS" (HH with Many Sodium 
slow inactivation processes, section 4.5.4). In the second case, 
those processes were added in a multiplicative manner (as dif- 
ferent processes affecting the same channel, in the uncoupled 
approximation), so s was replaced with J^[ ; s, in the voltage equa- 
tion (Equation 21). We denote this model as "Multiplicative 
HHMS" (section 4.5.5). In both cases, our analytical approxima- 
tions seemed to hold quite well. For example, the approximated 
Sy (/) (Equation 13) corresponded rather well with the numerical 
simulation of the full model (Figures 5B,D, respectively). 

Next, to test the limits of our assumptions we extended the 
HHS model to the HHSIP model (from Soudry and Meir, 2012b, 
see section 4.5.6) and added a potassium inactivation current 
which had faster kinetics (so r s ~ 5 Hz). So if Tt/ 1 = 10 Hz, we 
get T* «a 0.5r s , so the timescale separation assumption 2 is not 
strictly valid here. Also, for certain parameter values we get a limit 
cycle in the dynamics of s m , so the fixed point assumption 3 fails. 
However, it seems that our approximations still follow the numer- 
ical simulation of the full model: for p* at various stimulation 
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FIGURE 4 | Comparing mathematical results with full model simulation 

when the assumptions fail to hold. In the HHSIP model (HHS with 

potassium inactivation) we plot (A) p» (T^ 1 ) for different currents 

(to = 7.5, 7.7, 7.9, 8.1 , 8.3/iA from bottom to top). (B) Sy (f) for two values of 

7",. As before, "Sim" is a simulation of the full model, "Approx" is the 

analytical approximation, and " Ident" is the PSD Sy (f) of the linear system 

identified from the spiking data. Upper figure shows the case when 

T, « 0.5r s so the timescale separation assumption breaks down. In the lower 



figure the parameters are close to a Hopf bifurcation where a limit cycle is 
formed so the fixed point assumption breaks down, so the estimation of the 
limit cycle frequency component is less accurate. (C) The estimation of Si for 
T, -1 « 30 Hz is even better than in the HHS case. Similarly to (A-C) we plot 
the results of the HHMSIP model (HHSIP with many additional slow sodium 
inactivation kinetics) in (D-F), which has considerably more noise in the slow 
kinetics, and so even larger fluctuations (which further invalidates the fixed 
point assumption). See section 4.5 for various model details. 
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FIGURE 5 | Comparing mathematical results (green) with full model 
simulation (blue) for various models. (A) Coupled HHS (HHS coupled slow 
and rapid kinetics) (B) HHMS (HHS with many additional slow sodium 
inactivation kinetics) (C) HHSTM (HHS with a synapse) (D) Multiplicative 

frequencies X^ 1 and currents In (Figure 4A), for Sy (/) at T^ 1 = 
10 Hz when assumption 2 breaks down (Figure 4B, top), for 
Sy (/) at T^ 1 = 30 Hz when assumption 3 breaks down (near a 
Hopf bifurcation) and a limit cycle begins to form (see Figure 4B, 
bottom), and for state estimation of si using a linear optimal filter, 
again at T^ 1 = 10 Hz (Figure 4C). 

The only discrepancy seemed to appear in the limit cycle case, 
where the frequency of the limit cycle "sharpens" the peak in 
Sy (f) (Figure 4B, bottom). This may suggest that, in this case, 
the perturbations of the system near the limit cycle could be lin- 
earized, and that the eigenvalues of that linearized system might 
be related to the eigenvalues of the linearized system around the 
(now unstable) fixed point s*. More generally, the results so far 
indicate that even if our assumptions are inaccurate, it is possible 
that the resulting error will not accumulate and remain small - in 
comparison with the intrinsic noise in the model. 

Next, to challenge the approximation even more, we added to 
the HHSIP model four types of sodium currents with increas- 
ingly slower kinetics and fewer channels, similarly to the HHMS 
model (so this is the "HHMSIP" model, section 4.5.7). This 
significantly increased the variance of the dynamic noise n m , ren- 
dering the dynamics more "noisy." These random fluctuations 
in s m (Figure 4E) are of similar magnitude to the width of the 
threshold (non-saturated) region inpAP (s m ) (see Figure 6). This 
renders the fixed point assumption 3 inaccurate, since now the 
linear approximation pAP (s m ) ~ P* + w T s m breaks down most 
of the time. However, even in this case, the approximations seem 
to hold quite well with simulations of the full neuronal model 
(Figures 4D-F). 

In Figure 5A we used a coupled version of the HHS model 
("coupled HHS" model, section 4.5.2), in which the equations 
for r and s in the full model are tangled together, and not sep- 
arated as we assumed in Equations (2, 3). Even in this case, our 
approximations seemed to hold well. 

Finally, in Figure 5C, we extend the HHS model so that 
the stimulations are not given directly, but through a synapse. 



HHMS (variant of HHMS). As before, "Sim" is a simulation of the full model, 
"Approx" is the analytical approximation, and "Ident" is the PSD Sy (f) of the 
linear system identified from the spiking data. See section 4.5 for various 
model details. 

We used the biophysical Tsodyks-Markram model (Tsodyks and 
Markram, 1997) of a synapse with short-term depression, with 
added stochasticity ("HHSTM" model, section 4.5.3). This also 
seemed to work well. 

In all simulation we also added the PSD of the linear model 
identified from the spike output ("Ident."), to show that it can 
be estimated reasonably well. Note that the performance at the 
lowest frequencies seems to be significantly worse when they 
contain relatively high power. This is not surprising since it is 
typically harder to estimate model parameters, when the data has 
such (l// a ) PSD shape - which indicates long-term correlations 
(Beran, 1992). 

3. DISCUSSION 

In this work we found that under a temporally sparse ("spike- 
like") stimulation regime (Figures 1A,B) we can perform accurate 
semi-analytical linearization of the spiking input-output relation 
of a CBM (Figure 1C), while retaining biophysical interpretability 
of the parameters (e.g., Figure 7). This linearization considerably 
reduces model complexity and parameter degeneracy, and enables 
the use of standard analysis and estimation tools. Importantly, 
this method is rather general, since it can be applied to any 
stochastic CBM, with only a few assumptions. 

3.1. CONNECTION TO PREVIOUS WORK 

To the best of our knowledge, such results are novel, as no 
previous work examined analytically the response of general 
stochastic CBMs to temporally sparse input for extended dura- 
tions. However, the connection between sparse inputs and slow 
timescales has been previously made. It was previously suggested 
(Linaro et al., 2011) that sparse inputs could be used to identify 
neuronal parameters in a network of integrate and fire neurons 
with spike frequency adaptation. Interestingly, using different 
methods we reach a qualitatively similar conclusion here, though 
not in a network setting, and for a different class of neuron 
models. 
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Additionally, in Soudry and Meir (2012b) we modeled neurons 
under periodical stimulation using deterministic CBMs with slow 
kinetics, which are completely uncoupled from each other, and 
slower than the stimulation rate. Using a reduction scheme simi- 
lar in nature to that described here, we were able to describe the 
deterministic CBM's excitability and response using a discrete- 
time map - which "samples" the neuronal state at each stimula- 
tion. Analyzing this map, we obtained analytical results describing 
the neuronal activation modes, spike latency dynamics, mean fir- 
ing rate and short-time firing patterns. Stochastic CBMs were 
then examined numerically, and were shown to lead to qualita- 
tively different responses, which are more similar to the experi- 
mentally observed responses. 

The current work, therefore, generalizes this previous work. 
Here, we analyze the general case of stochastic CBMs, under gen- 
eral sparse stimulation patterns and with coupled slow kinetic 
dynamics. Therefore, the framework in the previous work 
(Soudry and Meir, 2012b) could be considered as a special case 
of this work, in which there is an infinite number of ion chan- 
nels (N — »- oo, so B r / S = D r / S = 0), T m = T* (so T m = 0) and 
A s (V) (the rate matrix) is a diagonal matrix. In the current work 
we similarly show that, in the generalized framework, the CBM's 
excitability and responses can be succinctly described using a 
discrete-time map. It is then straightforward to derive results par- 
alleling those in Soudry and Meir (2012b) in this more general 
setting, such as the mean firing rate (Equation 10). 

3.2. THEORETICAL NOVELTY 

However, the main novelty lies in our additional results, that 
could not be derived in Soudry and Meir (2012b). Specifically, 
due to the presence of noise, we were able to linearize the map's 
dynamics, and derive an explicit input-output relation. Such a 
linearization became possible because we made the (unusual) 
choice that the "input" to the CBM consists of the time-intervals 
between stimulation pulses, while the "output" is a binary series 
indicating whether or not an AP happened immediately after 
a stimulation pulse. The linearized input-output relation can 
be expressed either in biophysically interpretable "state space" 
(Equations 11, 12 and Figure 1C), or as a sum of the filtered 
input and filtered noise (Equation 20). Note that the overall I/O 
includes the mean output (Equation 10) which is non-linear. 
However, the linear part of the response, allows the derivation 
of the power spectral densities (Equation 13), the construc- 
tion of linear optimal estimators (e.g., Figure 3C) and blind 
identification of the (linearized) system parameters ("Ident." in 
Figures 3-5). 

Our results rely on three main assumptions. The temporal 
sparseness of the input tap <SC T m insures that the slow variables 
s m effectively represent the "neuronal state" alone (as V and r 
always relax to a steady state before the next stimulation is given). 
The additional assumption T m <^C r s allowed us to integrate the 
model dynamics and derive the reduced map (Equation 8) for 
the dynamics of s m , which is linear in T m . The last assump- 
tion is that the dynamics of s m can be linearized around a stable 
fixed point s*. This fixed point is generated due to the noisi- 
ness of the rapid variables (Figure 2), and the assumption T m <SC 
r s ensures that the stochastic fluctuations around s* are small. 



We performed extensive numerical simulations (section 2.5) that 
indicate that our analytical results are accurate - sometimes even 
if our assumptions break down. 

However, clearly there are cases, beyond our assumptions, in 
which are results cannot hold. For example, if T m has very large 
fluctuations, then the response of the neuron cannot be com- 
pletely linear, since 0 < Y m < 1. Such cases may require an exten- 
sion of the formalism described here. There are many possible 
extensions which we did not pursue here. For example, one can 
extend the modeling framework (e.g., multi-compartment neu- 
rons), stimulation regime (e.g., heterogeneous pulse amplitudes), 
or the type of neurons modeled (e.g., bursting and spontaneously 
firing neurons). However, it seems that an important assump- 
tion, that cannot be easily removed, is that the input is temporally 
sparse (r A p < T m ). 

3.3. PRACTICAL SIGNIFICANCE 

Is such a sparse temporally stimulation regime "physiologically 
relevant" for the soma of a neuron? Currently, such question 
cannot be answered directly, since it is impossible to accurately 
measure all the current arriving to the soma from the den- 
drites under completely physiological conditions. However, there 
is some indirect evidence. Recent studies have shown that the dis- 
tribution of synaptic efficacies in the cortex is log-normal (Song 
et al., 2005) - so a few synapses are very strong, while most 
are very weak. This indicates that the neuronal firing patterns 
might in fact be dominated by a small number of very strong 
synapses while the sum of the weak synapses sets the voltage 
baseline (Ikegaya et al., 2012). Such a possibility is supported 
by the fact that individual APs can trigger the complex network 
events in humans (Molnar et al., 2008; Komlosi et al., 2012). 
Also, in rats, individual cortical cells can elicit whisker move- 
ments in Brecht et al. (2004) and even modify the global brain 
state (Li et al., 2009). Taken together, these observations suggests 
that the above-threshold stimulation reaching the soma may be 
temporally sparse in some cases. 

There are other obvious cases were our results are immedi- 
ately applicable. First, in an axonal compartment, the relevant 
input current is indeed sparse - an AP spike train arriving from 
a previous compartment. Second, a direct pulse-like stimula- 
tion is used in cochlear implants (Goldwyn et al., 2012, and 
references therein). Lastly, such stimulation is used as an exper- 
imental probe (De Col et al., 2008; Gal et al, 2010; Gal and 
Marom, 2013). Specifically, since we now have a precise expres- 
sion for the power spectral density of the response, we are 
now ready to use these analytical results in Soudry and Meir 
(2014) to reproduce the experimentally observed l/f a behavior 
of the neuron and explore its implications on its input-output 
relation. 

4. METHODS 

In this section we provide the details of the results presented in the 
paper. Sections 4.1-4.5 here respectively correspond to Sections 
2.1-2.5. The first four (theoretical) sections can be read inde- 
pendently of each other (except when we discuss the repeating 
"HHS model" example). The last section give the details of the 
numerical simulations. 
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4.1. FULL MODEL (BIOPHYSICAL NEURON MODELS) 

As we explained in section 2.1, a general model for a biophysical 
isopotential neuron is given by the following equations 



V = f(V, r, 8 , /(*)), 
r = A r (V)r + B r (V,r)? r , 
s = A s (V)s + B s (V\ s)| s , 



(24) 
(25) 
(26) 



with voltage V, stimulation current I(t), rapid variables r 
(e.g., m, n, h in the Hodgkin-Huxley (HH) model Hodgkin and 
Huxley, 1952), slow variables s (e.g., slow sodium inactivation 
Chandler and Meves, 1970), rate matrices A r / S , white noise pro- 
cesses £ r / s (with zero mean and unit variance), and matrices B r / S 
which can be written explicitly using the rates and ion channel 
numbers (Orio and Soudry, 2012) (D = BB T is the diffusion 
matrix Gardiner, 2004; Orio and Soudry, 2012). In this section 
we give the specific forms of A r / S and B,-/ s , and their origin based 
on neuronal biophysics. 

Microscopic origins 

Such a model is commonly called a stochastic Conductance Based 
Model (CBM). In a non-stochastic CBM, the dynamics of the 
membrane voltage V (Equation 36) are deterministically deter- 
mined by some general function of V, the stimulation current 
1(f), and some internal state variables r and s. In contrast, 
the dynamical equations for r and s here adhere to a specific 
Stochastic Differential Equation (SDE) form, since these vari- 
ables describe the "population state" of all the ion channels in the 
neuron. We now explain the biophysical interpretation of those 
equations. 

At the microscopic level, each ion channel has several states, 
and it switches between those states with voltage dependent rates 
(Hille, 2001). This is usually modeled using a Markov model 
framework (Colquhoun and Hawkes, 1981). Formally, suppose 
we index by c the different types of channels, c = 1 , . . . , C. For 
each channel type c there exist channels, where each channel 
of type c possesses JC& internal states. In the Markov frame- 
work, for each ion channel that resides in state i, the probability 
that the channel will be in state j after an infinitesimal time dt is 
given by 



Af (V) dt 



A (c) (V) dt 



if j # i 
if j = i 



(27) 



where A^ (V) is called the "rate matrix" for that channel 
type. 

To facilitate mathematical analysis and efficient numerical 
simulation, we preferred to model stochastic CBMs using a com- 
pressed, SDE form. This method was initially suggested by Fox 
and Lu (1994), but their method suffered from several prob- 
lems (Goldwyn et al., 2011). In a recent paper (Orio and Soudry, 
2012) a more general method was derived, which had none of 
the previous problems, and was shown numerically to produce 
a very accurate approximation of the original Markov process 
description. 



Derivation 

According to Orio and Soudry (2012), if we define x^ to be the 
fraction of c-type channels in state k, and x' c) to be a column 



vector composed of Xu , then 



x (<0 = A M (y) x (c) + B (c) / V] x (c)\ £(c) 



(28) 



where f ^ is a white noise vector process - meaning it has zero 
mean and auto-covariance 



(f)(V c) (f')) T ) = I^5(f-f') 



where I is the identity matrix, S (f) is the Dirac delta func- 
tion, and S c ( j = 1 if c = d and 0 otherwise. Furthermore, B (c) 
is defined so that in Equation (28) each component of tf® , 
which is associated with a transition pair i ^ is multiplied by 

(a^x^ + Aj^Xij /N^ c \ and appears in the equation for x^ 

and & with opposite signs. Note that B^ is not necessarily 

square since it has rows but the number of columns is equal 
to the number of transition pairs. 

We now need to combine Equation (28) for all c to obtain 
Equations (1-3) . For simplicity, assume now that all ion 
channels types can be classified as either "rapid" or "slow" 
(this assumption can be relaxed). In this case we can concate- 
nate all vectors related to rapid channels r = y^\y ■ ■ ■ , x (j;)) 

and to slow channels s = fxjj +1 j, . . . , xjj +s ^ , where _R + 
S = C. We similarly define % r and % s together with the block 
matrices 



/A« 0 

0 A< 2 > 



A r 



0 0 



0 \ 
0 

A (R) ) 



/ A (R+D o 
0 A^+ 2 > 



\ 



V 



o 



o 



and similarly for B r and B s . Note that A r is square with 
size M = 5Zc= l rows while A s is square with size 
M=E?=/+ throws. 

4. 1. 1. Compressed formulation 

In some cases, it is more convenient to re-write Equations (1-3) 
in a compressed form (this is always possible) 



V = f(V, I, 8.1(f)), 

r = A r (V)r-b r (V) + B r (V, r)| f , 
s = A s (V) s - b s (V) + B 5 (V, s) § s , 



(29) 
(30) 
(31) 



where r, s, and £ r / 5 have been redefined (their dimension has 
decreased), as we will show next. First, we comment that the 
main disadvantage is of these equations is that they are less 
compact and the notation is somewhat more cumbersome. 
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However, there are also several advantages to this approach: 
(1) The vectors and matrices are smaller, (2) The rate and 
diffusion matrices do not have "troublesome" zero eigen- 
values and can be diagonal (which is analytically conve- 
nient), (3) Most CBMs are written using this form (e.g., the 
HH model), so it is easier to apply our results using this 
formalism. 

Derivation 

To derive these compressed equations, we use the fact x[ 

denote fractions, so XjfcXt. = 1> for all c. We can use this 
constraint, together with the irreducibility of the underlying 
ion channel process, to reduce by one the dimensionality of 
Equation (28) (see Soudry and Meir, 2012a for further details). 
Defining I to be the identity function, J to be the I with it 
last row removed, e = (0, 0, . . . , 1) T , u = (1, 1, . . . , 1) T , G = 
(I - eu T ) J T , A (c) = JA (C) G, B (c) = JB (c) (with x ( £ c) replaced by 

1 — x\ — %2 . . . — x K (c)_]) and = — JA^'e (A ° is invertible), 
we obtain the following equation for the reduced state vector 
y (c) = Jx (c) (which has only K {c) - 1 states) 



: A (c) y (c) -b + B W $M 



Again assuming that all channels can be classified as either 
'rapid" or "slow," we concatenate all vectors related to rapid 

,y^^ and to slow channels s = 
, y^+sj J , where R + S = C. We obtain Equations 



channels r : 



fa) 



(30, 31) by similarly defining b r ,b s ,£ r and £ s together with the 
block matrices 



/a (1) 



A,- 



0 A (2) 



0 0 



o \ 

0 



0 



A ( * +2) 



\ 



A (R+S) 



A s (V) = -y (V) - S (V) 
b s (V) = -S (V) 



(32) 
(33) 

B s (V,s) = V(S(V)(1 - s) + Y(V)s)Nr 1 4> (34) 
D s (V,s) = («(V)(1 - s) + y{V)s)N- 1 <t>. (35) 

Note that all the parameters are scalar in the HHS model, and so 
are not boldfaced, as in the general case. 

4.2. MODEL REDUCTION 

In this section we give additional technical details on section 
2.2. Specifically, we show how, given sparse spike stimulation 
and a few assumptions, it is possible to derive a simple reduced 
dynamical system (Equation 8) from the full equations of a gen- 
eral biophysical model for an isopotential neuron (Equations 
1-3), 



V=f(Y, r,s,I(f)), 
r = A r (V)r+B r (V, r)| r , 
s = A s (V) s + B 5 (V, s) $ s . 



(36) 
(37) 
(38) 



For more details on how its parameters and variables map to 
microscopic biophysical quantities, see section 4.1. 

4.2. 1. The excitability constraint 

As explained in section 2.1, we focus on models for excitable 
neurons describable by equations of the general form of 
Equations (36-38), rather than on arbitrary dynamical sys- 
tems. This imposes some constraints on the parameters 
(Soudry and Meir, 2012b). Formally, recall that tap and r s 
are the respective kinetic timescales of {V, r} and s, and 
that tap < t s . Suppose we "freeze" the dynamics of s (so 
that effectively t s =oo) and allow only V and r to evolve in 
time. We say the original model describes an excitable neu- 
ron, if the following conditions hold in this "semi-frozen" 
model: 



and similarly for B r and B s . Note that A r is square with 
M r = Ylc= l ^ <C ' — ^ rows while A s is square with M s = 
i — S rows. Furthermore, it can be shown (Soudry 

and Meir, 2012a) that A c is a strictly stable matrix (all its 
eigenvalues are also eigenvalues of A ( ^ except its zero eigen- 
value, and so have a strictly negative real part), and D c = 
B (l) B c is positive definite (so all its eigenvalues are real and 
strictly positive). Therefore, also A r and A s are both strictly sta- 
ble and D r and D s are positive definite. Therefore, if V is held 

constant, (s) and (r) tend to Soo = A s b 5 and = A r b,-, 
respectively. 

Example - the HHS model 

The HHS model can be easily written using the compressed for- 
mulation. For example, comparing Equation (23) with Equation 
(31) we find that 



1. If I (t) = 0, then for all initial conditions, V and r rapidly 
(within timescale tap) relax to a constant and unique steady 
state ("rest"). 

2. Assume that V and r are near rest, and a short stimula- 
tion pulse is given with duration f st ; m < tap and amplitude 
Jo- For certain initial conditions and values of in, we get 
either a stereotypical "strong" response ("AP response") or 
a stereotypical "weak" response in V ("no AP response"). 
Only for a very small set of initial conditions and val- 
ues of Iq, do we get an "intermediate" response ("weak AP 
response"). By "stereotypical" we mean that the shape of 
response does not change much between trials or for dif- 
ferent initial conditions in {V, r} (however, it can change 
with s). 

Note that due to condition 1, such an excitable neuron is not 
oscillatory and does not spontaneously fire APs. 
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4.2.2. Problem formulation 

Formally, suppose an excitable neuron receives a train of identical 
stimuli, so 



I(t) = £ n(t 



where n (#) is a pulse, of width t w (so n (x) =0 for x out- 



side [0, t w }). We denote by {F m }~ 
AP responses at times {f,„}^ = _ 00 
stimulation time t m (Figure 1A), 



_ 0O the occurrence events of 
i.e., immediately after each 



1 , if an AP occurs 
0 , otherwise 



(39) 



Defining T m = t m +i — t m , the inter-stimulus interval, and tap 
as the upper timescale of an AP event (Figure IB) we make the 
following assumption. 

Assumption 1. (a) The stimulation pulse width is small, t w < 
TAp.(b) The spike times {t m }^ =0 are temporally sparse, i.e., tap <3C 
T m for every m ("no collisions"). 

Our main objective here is to mathematically characterize 
the relation between {Y m } and {T m } under the most general 
conditions. 

4.2.3. Derivations 

We define the sampled quantities V m = V(t m ), r m = r(t m ), 
s m = s(f m ), x m = (y m ,r^,s^) T and the history set U m = 

{fofcJfc^oo ■ {Tk}k=-oo . ™ifc=-oc} ( n0te that H >n C K m + l). 

The Stochastic Differential Equation (SDE) description in 
Equations (36-38) implies that x m is a "state vector" with the 
Markov property, namely it is a sufficient statistic on the his- 
tory to determine the probability of generating an AP at each 
stimulation, 

P(Y m = l\x m ) = P(Y m = l\x m ,H m -i), (40) 
and, together with Y m and T m , its own dynamics 

P(x m +i|x m , T m , Y m ) = P(x m+ i|7Y m ) , (41) 
which implies the following causality relations 



To Ti 

Xfj -> Xl -> x 2 • • • x m 

I / I / I I ■ 

Y 0 Y x Y 2 ■ ■ ■ Y m 



(42) 



This causality structure is reminiscent of the well known Hidden 
Markov Model (Rabiner, 1989), except that in the present 
case the output Y m , affects the transition probability, and we 
have input T m . Theoretically, if we knew the distributions in 
Equations (40) and (41), as well as the initial condition -P(xo), 
we could integrate and find an exact probabilistic I/O relation 



p ({ y *}fc=u K^fcLo)- However, since it may be hard to find 
an expression for P(x m+ i|x m , T m , Y m ) in general, we make a 
simplifying assumption. 

Assumption 2. T m <g r s for every m. 

This assumption, together with Assumption 1 and the excitable 
nature of the CBM, renders the dynamics between stimula- 
tions relatively easy to understand. Specifically, between two 
consecutive stimulations, the fast variables (V (t) , r (t)) follow 
stereotypically either the "AP response" (Y m = 1) or the "no-AP 
response" (Y m = 0), then equilibrate rapidly (within time tap) 
to some quasi-stationary distribution q{V, r|s m ). Meanwhile, 
the slow variable s(f), starting from its initial condition at 
the time of the previous stimulation, changes slowly accord- 
ing to Equation (38), affected by the voltage trace of V (f) 
(through A S (V)). 

Summarizing this mathematically, we obtain the following 
approximations 



(43) 



P(Y m \s m )*> j P(Y m \V, r, s m )q(V, r \s m )dVdr, 

P ( + 1 1 + 1 j S/« + 1 1 Sm ■> Tm i Y m ) ~ <\ ( V m 

+ 1, fm+llSm+l) 

P(s m+1 |s m , T m ,Y m ). (44) 

Using these equations together with Equations (40) and (41), we 
obtain 

p A p(s m ) =P(Y m = l\s m ) = P(Y m = l\s m ,H m -i), (45) 
P (Sm + 1 \ *hn j ^ T m i T m ) = P (s m -|- i |s m , Y m , T m , T~Cm — l) • (46) 

Therefore, the "excitability" vector s m can now replace the full 
state vector x m = (V m , r^, s^) as the sufficient statistic that 
retains all relevant the information about the history of previ- 
ous stimuli. Given the input {T m }^ ) = _ 0O , Equations (45) and (46) 
together completely specify a Markov process with the causality 
structure 



To Ti 

So -> Si -> S 2 • • • S m 

I / I / I I 

Y 0 Yi Y 2 — Y m 



(47) 



Since the function pAP (s) is not affected by the kinetics of s, it 
can be found by numerical simulation of a single AP using only 
Equations (36, 37), when s is held constant (see section 4.2.4). 
Now, instead of finding P (s m+ i|s m , Y m , T m ) directly, we calcu- 
late the increments As m = s m + 1 — s m by integration of the SDE 
in Equation (38) between t m and tm+i. First, we integrate the 
"predictable" part of the increment 



(As m |s m , T m , Y n: 



L A 

a 



(V(t))s(t)dt, (48) 



A s (V (0) dt s„ 



(49) 
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to first order, where (X\Y) denotes the conditional expectation 
of X given Y. Note that A 5 ~ O (t^ 1 ), so second order correc- 
tions are of order O HX,^, -1 ) 2 } . Due to assumption 2, we have 
T m r s T 1 <SC 1, so these corrections are negligible. Now, 



/ A s (V (0) dt = r AP / A s (V (r)) dt 

Jt m V tap Jt m 



+ (T m -r AP )(— 1 f + ' A,(V(f))<ft 

\ 1 m — tap Jt m +TAp 

tap (A+ (s m ) Y m + A_ (s m ) (1 - Y m )) 
+ (T m - tap) A 0 (s m ) 



where we defined 
An (s m ) 

A— (s m ) 



1 rhn + i 

/ A s (V(t))dt 

Jt m + XkV 



T m — TAP Jf m +r A p 

TAP A m 
1 |-<m+rAP 



(50) 

ify m = o, (5i) 



A+(s m )=— / A 5 (V(t))df, ify m =l, (52) 

tap Jt m 



which are the average rates during rest, during an AP response 
and during a no-AP response, receptively. Note a similar notation 
was also used in Soudry and Meir (2012b) (e.g., Equations 2.15, 
2.16 there), where the +/ — /0 were replaced with H/M/L. 

Next, we calculate the remaining part of the increment, which 
is the "innovation," 

ti m = As m (As m |s m , T m , Y m ) . 

Obviously, (n m |s m , T m , Y m ) = 0, and also 



(n m n |s m , T m , Y n 



^ m+1 B s (V(f),s(f))| 5 (f)^ 



atm+l -,- 
B s (V(f),s(f')) t- s (t')dt'y \s m ,T m ,Y„ 

j m + 1 dtj^ dt'S(t - t')B s (V (f) , s(f)) Bj (V (t') , s (f')) 

J m + 1 B s (V (f) , s (0) Bj (V (f') , s (0) 

^'" + 1 D 5 (V(f),s(f))df 

f tm + 1 

L D - 



; (V(t),s m )dt 



to first order. Note that D s ~ O (r s V w )» where N = min c N (c) 
(N^ is the channel number of the c-type channel, as we defined 
in section 4.1), while Equation (53) has corrections of size 
0^(r m T"Vw) 2 ). Since N> 1 (usually N » 1), and due to 



assumption 2, we have T m r s /N <C 1, so these corrections are 
also negligible. Now, 



L 



D 5 (V (f) , s m ) dt 



(53) 



tap (Y m D + (s m ) + (1 - Y m ) D_ (s m )) + (T m - t A p) D 0 (s m ) 



where we defined 
Do (s m ) 

D- (Sm) 



1 ftm + 1 

/ D s (F(f),s m )dr 



T m T A p Jt m +Tjtf 
tap Jt m 



(54) 



(y(r),s m )df, if 7 m = 0 (55) 



/•Im + IAP 

D+ (s m ) = D 5 (V (0 , s m ) dt, if Y m = 1. (56) 

tap Jt m 

Additionally, we note that A±/o (s,„) generally tend to be rather 
insensitive to changes in s m . This is because the kinetic transi- 
tion rates (which are used to construct A s (V), as explained in 
section 4.1) tend to demonstrate this insensitivity when simi- 
larly averaged (see Figures 4B, 5 in Soudry and Meir, 2012b). 
The usual reasons behind this are (see appendix section Bl of 
Soudry and Meir, 2012b): (1) The common sigmoidal shape of 
the voltage dependency of the kinetic rates reduces their sensitiv- 
ity to changes in the amplitude of the AP or the resting potential 
(2) The shape of the AP is relatively insensitive to s (3) The 
resting voltage is relatively insensitive to s. Therefore, in most 
cases we can approximate A±/o (s m ) to be constant for simplic- 
ity (though this not critical to our subsequent results), as we shall 
henceforth do. 

Additionally, we note that, strictly speaking, the voltage trace 
during an AP and at rest are stochastic, and therefore, A+, A_, 
An, D + , D_ and Do are stochastic. However, there are two fac- 
tors that render this stochasticity negligible. First, the sigmoidal 
shape of the kinetic rates implies that A(Y^) is rather insensitive 
to fluctuations in the voltage (Figure 4 in our Soudry and Meir, 
2012b). Second, noise mainly plays a role in the timing of AP 
initiation, but does not much affect the AP shape above thresh- 
old (see AP voltage traces in Schneidman et al, 1998, p. 1687). 
Therefore, we shall approximate A + , A_, Ao, D + , D_ and Do to 
be deterministic. 

In summary, defining 

A(Y m , T m ) = r AP T m l (Y m A+ + (1 - Y m ) A_) + 

(1 - tapT^Ao, (57) 



and 



D (Y m , T m , s m ) = T A pT m (Y m D + (s m ) + (1 - Y m ) D (s m )) 



we can write 



-^(I-tap^^Do (s m ) 



As m — T m A. {Y m , T m ) s m + n m , 



(58) 



(59) 
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with(n m |s m , T m , Y m ) = 0 and 

/n m n m |s m , T m , Y m ^j = T m D (Y m , T m , s m ) . (60) 

These equations correspond to the result presented in 
Equation (8). 

Finally, we note that the distribution of n m given s m , T m , Y m 
can be generally computed using the approach described in Orio 
and Soudry (2012). For example, it can be well approximated to 
have a normal distribution if channel numbers are sufficiently 
high and channel kinetics are not too slow (Orio and Soudry, 
2012). In that case only knowledge of the variance (Equation 60) 
is sufficient to generate n m . And so, using Equations (45), (59) 
and the full distribution of n m , we can now simulate the neuronal 
response using a reduced model, more efficiently and concisely 
(with fewer parameters) than the full model (Equations 36, 38), 
since every time step is a stimulation event. The simulation time 
should shorten approximately by a factor of (T m ) /dt, where dt 
is the full model simulation step. Note that the reduced model 
parameters, having been deduced from the full model itself, still 
retain a biophysical interpretation. 



4.2.4. Calculation o/p AP (s) 

We numerically calculated pAP (s) by disabling all the slow kinet- 
ics in the model - i.e., we only use Equations (1, 2) in main text, 
while s = 0. Then, for every value of s we simulated this "semi- 
frozen" model numerically by first allowing r to relax to a steady 
state and then giving a stimulation pulse with amplitude Iq. We 
repeat this procedures 200 times for each s, and calculate pAP (s) 
as the fraction of simulations that produced an AP. A few com- 
ments are in order: (1) In some cases (e.g., the HHMS model) 
we can use a shortcut and calculate Pap (s) based on previous 
results. For example, suppose we know the probability function 
Pap (s) for some model with a scalar s and we make the substitu- 
tion s = h (s) where the components of s represent independent 
and uncoupled channel types (Orio and Soudry, 2012) - then 
Pap ( s ) = Pap (h (s)) in the new model. (2) The timescale sepa- 
ration assumption tap <3C T m r s implies that all the properties 
of the generated AP (amplitude, latency etc.) maintain similar 
causality relations with s m as does Y m , so we can find their dis- 
tribution using the same simulation we used to find Pap (s), 
similarly to the approach taken to compute L (s) in the deter- 
ministic setting (Soudry and Meir, 2012b). (3) Numerical results 
(Figure 6) suggest that we can generally write 




0.95 



0.85 

4.6 
4.5 
4.4 









"*°-. 


""©.._ 




xKT 3 


o 






o 


0 


o 

o 



8.2 8.4 



k)g 10 JV = 4 
log 10 AT = 5 
N = 6 
AT = 8 

log^A^lO 

| login N =12 

0.94 0.96 




10" 
10" 

i 

10" 
10" 



o. 


■"hh 


fN fit 


•a. 




















**©., 








""O 



10 10 



10 10 

N 



10 10 



FIGURE 6 | Fitting of p flP (s) = 4> ((s— a) /b) in the HHS model. (A) Fitting of pap (s) for various values of Iq. (B) Fitting shows that a is linearly decreasing 
in Iq. (C) Fitting of p AP (s) for various values of N. (D) Fitting shows that b oc 1/VW. 
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FIGURE 7 | The averaged kinetic rates. Left: The averaged rates 
demonstrated for three common kinetic rates y (V) with sigmoidal shapes. 
Right: The voltage threshold of the sigmoid determines whether the process 
is sensitive to APs (the output), stimulation pulse (the input), or neither. Note 
that a similar classification of biophysical processes affecting excitability was 
previously suggested inWallach (2012, Figure 3.1 ). 



p A p(s) = 4>(e(s)/Jn?), (61) 

where O is the cumulative distribution function of the nor- 
mal distribution, E (s) is some "excitability function" (as defined 
in Soudry and Meir (2012b), so pAP (s) = 0.5 on the threshold 

— 1/2 « , 

0 = {s\E (s) = 0}), and N r , the "noisiness" of the rapid sub- 
system, directly affects the slope of Pap (s) (Figure 6D, bottom). 
Also, as explained in Soudry and Meir (2012b), £(s) is usually 
monotonic in each component separately and increasing in Iq 
(Figure 6C, top) - which could be considered as just another 
component of s which has zero rates. 

4.2.5. Compressed formulation - reduction 

We can perform a very similar model reduction and lineariza- 
tion using the compressed formalism presented in section 4.1.1. 
We just need to define (or re-define) A±_n, b±,o> O±,o (s m )> 
A (Y m , T m ), b (Y m , T m ) and D (Y m , T m , s m ) in the obvious way 
and repeat very similar derivations, arriving to 

As m = T m ^A (Y m , T m ) S m b (Y m , ^m)] ~t~ rim> 

instead of Equation (59) (or Equation 8). Next, we demonstrate 
this for the HHS model. 

4.2.6. Example - HHS model reduction 

We derive the parameters of the HHS reduced map. Recall that the 
HHS model is based on the compressed formulation. Following 
the reduction technique described in the previous sections, we 
numerically find the average rates y±,o an d <5± o (as in Equations 
(2.15, 2.16) of Soudry and Meir (2012b), where there we denoted 
H/M/L instead of +/ — /0 here), tap andpAp (s) (section 4.2.4). 
From Equations (32, 35), we find, 

■A±,o = -y±.o - 5±,o (62) 
b±,o = S±,o (63) 
At.o 0) = K l (a ±i0 (1 - s) + y±,qs) . (64) 

and so A ( Y m , T m ) and D(Y m , T m , s m ) are defined as in Equations 
(57) and (58), and similarly 

b(Y m , T m ) = r AP T m 1 (Y m b + + (1 - Y m ) bJ) +(l - Tap J" 1 ) b 0 . 

We give for example some specific values: if tap = 15 ms, then 
in the range Iq = 7.5 — 8.3 fiA, we have 5±,o = 25.5 — 25.6 mHz, 
Y+ = 22.9 - 22.1 mHz, y_ = 0.9 - 1.3/x Hz and y 0 = 0.29 - 
0.28/x Hz. 

Recall that these averaged kinetic rates are determined by 
the shape of the voltage dependent rates (y (V) and S (V), see 
Equation 125) (Soudry and Meir, 2012b). The relative values of 
the averaged kinetic rates determine what kind of information 
can be stored in s (which retains the "memory" of the neu- 
ron between stimulation). We qualitatively demonstrate this in 
Figure 7 depicting the values of y±,o for three different shapes 
of y (V): when y (V) is sigmoidal with high threshold, when it is 
sigmoidal with low threshold and when it is constant. These deter- 
mine whether y (V) is affected by the output (APs), the input 



(stimulation pulses) or neither. Therefore: (1) if y (V) and S (V) 
are independent of the voltage, then s cannot store any informa- 
tion on input or output. (2) if y (V) or S (V) have low voltage 
threshold , then s can directly store information on the input. (3) 
if y (V) or S (V) have high voltage threshold , then s can directly 
store information about the output. In the HHS model the inacti- 
vation rate y has high threshold, while S is voltage independent - 
therefore, s directly stores information on the output. 

4.3. LINEARIZATION 

In this section we present a more detailed account on how 
to arrive from the reduced model (mainly, Equation 8) to its 
linearized version (the results in Equations 11, 12). 

First, we write the complete reduced model, using Equations 
(59), (60), and (45). The reduced model is a non-linear stochastic 
dynamic "state-space" system with T m , the inter-stimulus interval 
lengths, serving as inputs, s m representing the neuronal state, and 
Y m the output. We have 

Asm — T m A (Y m , T m ) Sffj -\- n m , (65) 
Y m = Pap (s m ) + e m , (66) 

where ^n m n m |s m , T m , Y m ^j = T m D (7 m , T m , s m ), 

A(Y m , T m ) = r A? T m l (Y m A+ + (1 - Y m ) A_) + 
(I-tapT-^Ao, 

and 
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D(Y m ) = r AP T m ' (Y m D + (s m ) + (1 - Y m ) D (s m )) + 



and we defined 



(l - tapT^Do (s m ). 



— Y m Pap (fhn) • 



(67) 



Based on the causality structure in Equation (42), it is straight- 
forward to prove that e m and n m are uncorrelated white noise 
processes - i.e., (e m ) = 0, (n m ) = (e„n m ) = 0 and (n,„n^ = 
(n m nT)3 m „, (e m e„) = {e 2 m )S mn where S mn = 1 if n = m and 0 
otherwise. 

We now examine the case where {T m } is a Wide Sense 
Stationary (WSS) process (i.e., the first and second order statis- 
tics of the process are invariant to time shifts), with mean T*, 
so that the assumptions tap « T m « t s are fulfilled with high 
probability. In this case the processes {s m } and {Y m } are also 
WSS, with constant means (s m ) = s* and (Y m ) = p*. Also, it is 
straightforward to verify that /r m n„j = 0, and |r m e„j = 0. 

In order to linearize the system in Equations (59-66) 
we denote T m = T m — T*, Y m = Y m — p*, s m = s m — s*, w = 
VpAp| s »- In order for this linearization to be accurate we require 
that s„, is "small enough." 



Assumption 3. 

wise) and \ w T s t 



With high probability 
, (VVpAp| s »)s« 



» si 



Sm \ <3£ I s * I (component- 



This assumption essentially means that s* = s* (p*, T*) is a stable 
fixed point of the system (Equations 59-66), and stochastic fluc- 
tuations around it are small, compared to the size of the region 
{s|pAP (s m ) 7^ 0, l} (usually determined by the noise level of 
the rapid system {V, r}, see section 4.2 A). Note that the region 
is usually rather narrow (Figure 6) and therefore |s m <SC |s*| is 
often implied by this description. Given Assumption 3, we can 
approximate to first order 



PAP (Sm) 



: p* + W T S„ 



(68) 



which allows us to linearize Equation (66). This essentially means 
that the components of s m determine the neuronal response lin- 
early, with the components of w serving as the effective weights 
(related to the relevant conductances in the original full neuron 
model). 

Next, we wish to linearize Equation (59). Using our assump- 
tions, we obtain to first order 

Sm+i ^ s m + A(p*, T*) (s* +s m ) (69) 
+ A 0 (s» + s,„) f m + tap (A+ - A_) (s, + s m ) Y m + n„, 
Taking expectations and using Equations (66) and (68), we obtain 

o = (s m+ i -s m ) » r*A(p*, r»)s*, (70) 

to zeroth order. Defining the solution of this equation is 
s* (p*, X*) and we can findp* implicitly from 



We write the explicit solution of this equation as p* Next, 
using |s m | <C |s*|, Equation (70) and defining 



F = I + T* A* (p* , r„) 
d = Aos* 

a = tap (A + - A_) s* 
we can approximate Equation (69) as 

l = Es m -f- dT m + &Y m + n n 
which, together with 



Y m — w s m + e m , 



(72) 
(73) 
(74) 



(75) 



(76) 



yields a simple linear state space representation with T m as the 
input, s m as the state, Y m as the output and two uncorrelated 
white noise sources with variances 



to first order. 

4.3. 1. Derivation of w 

From Equation (61), we note that generally we can write 



(77) 
(78) 



w = VpAP (s) s 



exp ( ) , (79) 



2N r 



where in many cases the excitability function E (s) has the 
form E (s) = /t T s — 8, where the components of fi are propor- 
tional to the relevant conductances (Soudry and Meir, 2012b). 
Therefore, if 

P* = Pap (s*) = <$ (E(s) Oor 1 

then E (s) — ► ±oo, so in this case (assuming E (s*) is not a 
particularly "pathological" function) we have 



0. 



(80) 



4.3.2. Compressed formulation - linearization 

In the compressed formulation (introduced in sections 
4.1.1 and 4.2.5), we can perform similar linearization by 
re-defining F = I + T*A (p*, T*) , d = A 0 s* - b 0 , a — tap 
((A + — A_) s* — (b+ — b_)), and repeat very similar 
derivations, where now we can write more explicitly 



s* = A 1 (p*, T*) b (p*, T*) . 



(81) 



P* = Pap (s* (p*, 



(71) instead of Equation ( 70 ) . 
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4.3.3. Example - HHS model linearization 

Note again that all the parameters are scalar now, and so are not 
boldfaced, as in the general case. From Equations (71) and (81) we 
obtain s* and for a given T* . Once s* is known, from Equation 
(79) w can be obtained 7 . Next, we denote the average inactivation 
rate at steady state by 



Y* 



{p*Y+ + (l - P*) Y-) tApT* 1 + (l - r A pi; ') Yo, 



and similarly for the recovery rate 5*. And so, s* = <5*/ (y* +5*), 
and 

A* = A* (p*, s*) = — y» - 5*, (82) 
fo* = -5*, (83) 
-D* = D* (p„, T», 5*) = AT 1 (5,y,/ (y, + 5*)) . (84) 

Denoting yi = y+ — y_ and similarly for Si, we obtain 

J? = 1 - T* {y* + &*) (85) 
« = tap (X*<5i - yi<S*) / (y* + 5*) (86) 
d = (Y*$o ~ Yo&*) I (Y* + 8*) (87) 
Finally, from Equations (77, 78) we find 

E„ = T*D* (88) 
°e=P*-pl (89) 
4.4. LINEAR SYSTEMS ANALYSIS 

In section 2.4 we describe the neuronal dynamics using a linear 
system for the fluctuations, as depicted in Figure 1. This linear 
description allows us to use standard engineering tools to analyze 
the system. In this section we provide an easy to follow description 
on how this was done, for those unfamiliar with these topics. 

4.4. 1. Second order statistics and linear systems 

We start with a short reminder on some known results for 
stochastic processes (Papoulis and Pillai, 1965; Gardiner, 2004); 
these results are standard but are provided for completeness. 
These results will be used in later sections. 

Assume {x„,} and {y m } are two real-valued vector stochastic 
processes that are jointly wide-sense stationary (i.e., a simultane- 
ous time shift of both processes will not change their first and 
second order statistics). We define the cross-covariance (recall 
that x = x — (x) ) 



and the Cross-Power Spectral Density (CPSD), given by its 
Fourier transform 



S xy (co)±F[R 1SY (-)](co) = 



—io)k 



k=—oo 



7 Also, as explained in section 4.3.1, we approximately have w oc gNa> from 
Equation (21). 



Additionally, the auto-covariance is defined as R x = R^ and the 
corresponding Power Spectral Density (PSD) as S x = S^. Also, 
note that ,Ryx (k) = R^y i—k) and so Syx (k) = S^y (—<«). 

Suppose now that {y m } is generated from a process {x m } 
using a linear system: i.e., if the Fourier transform x (w) = 
X!|fcl-oo %ke~ exists, then in the frequency domain 

y {w) = H (o>) x (w) , 

where H (&>) is a matrix-valued "transfer" function. Therefore, 
under some regularity conditions (allowing us to switch the order 
of integration end expectation), 



Sxy (co) 



k = — oo 



And similarly 



S Y (w) 



S X (»)H' (w) 



OO 

k= — oo 



(90) 



-icok 



H (—ft)) S x (ft>) H (<m) 



(91) 



where in the second equality here we used an almost identical 
derivation as for Sxy (ft)). 
Note that if instead 

y (co) = H x (co) x (ft>) + H z (ft>) z (to) , 

where x and z are two uncorrelated signals, then we can write 

y (ft)) = H (ft)) v (cd) , 

where 



H (ft>) 



H x (co) 0 
0 H z (ft)) 



V(ft)) = [X(ft)),z(ft))] . 



Thus Equations (90) and (91), respectively give 

Sxy (CO) = S x (CO) Hj (ft)) , (92) 

S y (ft)) = H x (-«) S x («) Hj (co) + H z (-«) S z (co) Hj (&>). (93) 

4.4.2. Ffte second order statistics of our system 

Previously, we derived Equations (11, 12), which describe the 
neuronal dynamics using a linear system, written in "state-space" 
form 



Sm+l 



Fs„, + dT m + aY m + n„ 

w T s,„ + e m 



(94) 
(95) 



where n m , e m and T m are uncorrelated, zero mean processes with 
the PSDs E n = T^D (p*, T*, s*) , of = (l - p*) and Sj («), 
respectively. 
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In order to apply Equations (92) and (93) to our system we Note that if the dimension of s is finite and there is no 
first need to find the transfer function of the system. Applying the degeneracy, we can always write 
Fourier transform to Equations (94, 95) gives 

M _ 

e !a, s(&>) = Fs(a>)+dr(«) + aF(ft>)+n(ft>), (96) S F (f) = c 0 + ^] ' (107) 

Y (to) = w T s (to) + e (to) . (97) 

where X;, the poles of Sy (/), are determined solely by the poles 
Re-arranging terms, we obtain of Hc (y) and Sj . ^ while a u the other parameters in Equation 

(106) affect only the constants Cj. Commonly, Sy if) has no 
s (ft>) = H c (w) (ft)) + df (ft)) + ae (to)) , (98) poles - for example, if Sj if) is constant so T m is a renewal process 

(e.g., the stimulation is periodic or Poisson). Therefore all poles 
Y (a,) = w T H c (ft.) (n(ft>) + df (ft)) + ae (ft))) + e (to) , (99) of Sy (/) (or the other PSDs) are determined by H c if), i.e., Xj are 

the roots of the characteristic polynomial 

where we denoted , 

AI — A (p* , X*) — X^ 1 aw T |=0. (108) 

H c (co) = (le ,ffl -F-aw T ) _1 . 

4.4.3. Spectral factorization 

This gives the "closed loop" transfer functions of the system Equations (96) and (97) can be re-arranged as a direct I/O rela- 
(including the effect of the feedback Y(co)). Next, combining tion > formulated, for convenience, in the frequency domain (this 
Equations (98, 99) and Equations (92, 93 ) leads to explicit expres- can be either / or w ~ in the section we use w for brevu 7 of nota " 
sions for the PSDs and CPSDs ^on, an< ^ / m omer places). Specifically, this relation is of the 

form 

S sT («) = H c (-ft)) dS T (ft)) (100) „ ( 

Y(fij) = H ext (ft))T(ft))+H mt (ft))v(ft)), (109) 

S s M = H c (-w)(T n + aa T a 2 + dd T SrMW (to), (101) 

so v m = T 1 (v (ft))) is a single scalar "noise" process with zero 
Syr (ft>) = w T H c (—to) dSj (to) , (102) mean and PSD tr 2 (here T~ l is the inverse Fourier transform). 

/ T \ This v m process combines the contributions of e m and n m , which 

Sy(to) = w H c (-to) \Y, n + dd S r (to)jH c (to)w (103) are me noise p rocesse s in the original system (in Equations 

, 96, 97). Such a description, as in Equation (109), describes con- 

9 I T 

+ tr £ : 1 + w H c (—ft)) a . cisely the contributions of the input and noise to the output (an 

ARMAx model Lejung, 1999). Using 92 and 93 we respectively 

For low frequencies it is sometimes more convenient to nn<1 mat 
use the "continuous-time" versions of the PSDs, Sxy (f) = 

T.Sxy (to) w=27r/Tt for/ « T-\ which are approximated by Syt ^ = ^ Sr ^ (110) 

Ssr if) = T-H C (-/) dS r (f) * ^ = ^ (W) | 2 Sr (W) + | H ' nt H' * (1U) 

S s (f) = H c (-/) (D (p. , X* , s*) + 77 1 aa T tr 2 + X~ 2 dd T Comparing Equation (102) with ( 1 10) we obtain 

S T if))uJif), (104) H ext (to) =w T H c (to) d. (112) 

Syt if) = 27V 1 " H c (-/) dS r (/) , (105) Comparing Equation (103) with (111), while using Equation 

(112), will yield the equation 

s Y if) = w t h c (-/) (d (p„ n, «,) + r- 2 dd T s T if)) 

H J(f) w (106) |H int (to)| 2 a v 2 =w T H c (-to) EnHjMw + a 2 |1 

+ r.tr 2 |l + 77 1 w T H c (-/) a!' . + wT H C (-*>) a f • (H 3 ) 



wnere This is a "spectral factorization" problem (Anderson and Moore, 

1979), with solution 

H c (/) = (2^I-A(p„r,)-r- 1 aw T ) \ H' nt M = w T H c (to) K+l, (114) 

and we used the fact that F = I + X* A (p* , X*) (Equation 72) and where 

E n = T*D(p*, T*,s*) (Equation 77). K = a + FPwtr v 2 (115) 
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and yield a method of solving the first problem (e.g., section 3.3 in 

a 2 = w T Pw + a 2 , (116) Ande rson and Moore, 1979). 

A relatively simple way to approach the second (filtering) 
with P the solution of problem involves the output decomposition we have found in 

section 4.4.3 

P = FPF T - (w T Pw + a 2 ) FPww T PF T + E n , (117) „ T / _ \ 

V / Y(co) = w 1 H c (co)dT(co) + (w T H c ((a)K + lj v(co) . 

derived from the general discrete-time algebraic Riccati equation. 

This can be verified by substitution ' Usin g this decomposition we can now write a new state-space 



representation for the system in terms of new state variable z„ 

w T H c (-«)£„H7Mw+ff e 2 |l+w T H f («)ar , s 

I I z m+ i = ^F + aw J z m + dT m + Kv,„, 



|H mt (co) | a 2 y m = w T z m + v m , 



(l) 



|2 

w T H 0 (co) FPwcr^ 2 + 1 a 2 



(2) 



n I — i j — 

and the estimation error is simply 



w T z„ 



w T H c (-co) - FPF T + <7,7 2 FPww T PF T ^ Hj (co) w which has the same output in the frequency domain (recall, from 

linear systems theory, that a single I/O relation can be generated 

+ a 2 i _|_ W T H C (w) a by multiple state space realizations). This "innovation form" is 

particularly useful, since, given the entire history of the previous 

It I 2 a f ~ ~ 1 m 1 

w H c (ft>) (a + FPwaC 2 ) + 1 cr 2 inputs and outputs H m _ i = I Tj., % \ , we can recursively 

I I J k=— oo 

estimate the current state precisely (with zero error) (Anderson 

[w T H„ (-ft)) (P - FPF T + cr i r 2 FPww T PF T J Hj (o) w + cr 2 and Moore, 1979) 

1 -w T H 0 (ft))a| z m = ^F + aw T ^z„,_i + dT m _i +K^7 m _i -w T z m _^ . 

[w T H 0 (—ft)) (P — FPF yN ) H T (ft>) w — w T Pw Given this precise estimate of z m , the best linear estimate of Y m is 

^ ' simply 

w T H 0 (w) FPw - w T H 0 (-«) FPwj 1 1 - w T H 0 (w) a| 2 

[w T (FH 0 (ft,) + I) P (F T Hj (a,) + 1^) w + w T H 0 (-w) FPF 
(ft)) w — w 1 Pw 

- w T H 0 (ft)) FPw - w T H 0 (-ft)) FPwj 1 1 - w T H 0 (co) ap' \( Ym ~ [ Y ^ H ^)) ) = ( V '«> = ff v ■ 

= " Since both the innovation form and the original form have the 

same second order statistics for the input-output, the optimal lin- 

where in (1) we used^ the fact that w T H c (ft)) = ear estimator (and its error) for Y m in the original system would 

w T H 0 (ft))(l -w T H 0 (ft))a) from the Sherman-Morrison be the same. Moreover, one can show (Anderson and Moore, 

lemma, with H 0 (w) = (e m l - F) _1 being the "open loop" 1979) that Equation (118) will also give the optimal linear esti- 

version of H c (co) (i.e., if a was zero), and in (2) we used the fact mate of s m in the original system, and with error P (Equation 

that H 0 (ft)) = e~ m (FH 0 (co) + I). 117). This solution is the well-known "Kalman filter." 

4.4.4. Optimal linear estimation of linear systems 4AS - Example - HHS model power spectral densities 

Given that the neuronal dynamics are given by the linear sys- Substituting the parameters for the linearized map (Equations 

tern in Equations (96, 97), there are two different estimation 85 " 89 ) int0 the expressions for the power-spectral densities 

problems one may be interested in. We may want to esti- (Equations 104-106), gives 
mate, based on the history of the previous inputs and outputs 

\f k ,Y k } m l .either the parameters ofthe model (F, w, a, d, a e r /rS w2 ( D * + ^^^(f^ + ^CT 2 ^/) 2 + A 2 ) 

I J*:=-oo JF(/) = j (11") 

and S n ), or the variables in the model (Y m or s m ). The first (2^/) 2 + (a* + T* 1 wa^J 
problem is generally termed a "system identification" problem 

(Lejung, 1999), while the second is a "filtering" (or predic- , , _ D* + T^ 1 a 2 a 2 + T~ 2 d 2 S T (f) 

tion) problem (Anderson and Moore, 1979). Both are intimately 5 V)~ . , ^ >2 

related, and sometimes the solution of the second problem can \ n t) \ * * wa J 
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S YT {f) 



Tl l wd 



lit ft ■ 



T± wa 



St (J). 



(121) 



Note that when St (f ) = 0 (i.e., periodical spike stimulus), Sy (f 
has the shape of high pass filter (Figure 3B, top). In contrast, S s (/ 
(Figure 3B, bottom) and Syr if) both have the shape of a low 
pass filter (Figure 3D, top). From Equations (111) and (110) we 
know that S Y (f) = \H int (f)\ 2 a 2 v and S YT if) = H ext (f) S T if), 
respectively. Therefore, this indicates that H mt (/) and H ext (/) are 
high pass and low pass filters, respectively. 

4.4.6. Power spectral densities of response features 

So far we have concentrated on the PSD of the response Y m . 
However, it is easy to extend our formalism to derive the PSDs 
of different features of the AP, such as its latency or amplitude. 
We exemplify this on the latency. In Soudry and Meir (2012b) we 
showed (Figure 3) that for deterministic CBMs, the latency of the 
AP generated in response to the m-th stimulation can be writ- 
ten as a function of the excitability L m = L (s m ). In a stochastic 
model, we have instead 



L (s m ) + <t> m 
not defined 



Ym 
Y 

> 1 n. 



where <p m is a zero mean, white noise process generated by the 
stochasticity of the rapid system. Since it is problematic to define 
the PSD of L m if sometimes Y m = 0, we focus on the case that 
pa, = 1, so we always have Y m = 1 . In this case, assuming again 
that the perturbations in s m are small, we can linearize 

L (Sm) ^ L (s*) + 1 s m 

where 1 = VL (s) s=Sj , to obtain (using Equation 11) 



2012b). The results are not affected qualitatively by our choice 
of a square pulse shape. We define an AP to have occurred if, 
after the stimulation pulse was given, the measured voltage has 
crossed some threshold (we use V t ^ = — lOmV in all cases). 
In all cases where direct stimulation is given, unless stated oth- 
erwise, we used periodic stimulation with In = 7.9 fiA and T* = 
50 ms. Note that for the parameter values used, no APs are spon- 
taneously generated, consistently with experimental results (Gal 
etal, 2010). 

The PSDs were estimated using the Welch method and aver- 
aged over eight windows, unless 1 // behavior was observed, in 
which case we used a single window instead, since long term cor- 
relations may generate bias if averaging is used (Beran, 1994). 
Numerical estimation of the cross-PSD is more problematic. 
When estimating cross-spectra, estimation noise level can be 
quite high (proportional to the inverse coherence, according to 
Bendat and Piersol, 2000, p. 321). To estimate the level of esti- 
mation noise, we estimate the cross-spectrum with the input 
randomly shuffled (Figure 8). Since in this case there is no input- 
output correlation, this new estimate is pure noise. Finally, as 
suggested by the reviewer, we smoothed the resulting PSD (or 
cross-PSD) in all figures (except in Figure 8, where we aimed to 
show the level of estimation noise). To achieve uniform variance 
with low bias, we divided the spectrum into 30 logarithmically 
spaced segments from / = 10~ 3 Hz to the maximal frequency 
{Tit/2). In each segment n the PSD (or cross-PSD) was smoothed 
using a window of size n. 

Next, we describe the models used Figures 3-5 and provide 
their parameter values. These models have either been studied 
in the literature or are extensions of such models, which are 
meant to explore the limit for the validity of our analytic approx- 
imations. In all cases where direct stimulation is given, unless 
stated otherwise, we use periodic stimulation with Iq = 7.9 [ih 



s m + 1 



Fs^ ; -\- dT m -j- n n 

1 Sm ~r~ 4>m 



(122) 
(123) 



where he F = I + T^A (1, T*). Therefore, it is straightforward to 
show that the PSD of the latency is 



S L (f)=l T S s (f)\+T^ 



(124) 



where er^ = [4> m ). Note that if latency is a good indicator of 
excitability, i.e., L (s) changes similarly to p (s) so that 1 oc w, then 
Si (/) = ciSy (f) + C2 for some constants c\, ci, when the input 
is periodic (T m = T*) andp* — >■ 1. 

4.5. NUMERICAL TESTS 

MATLAB (2010b) code is available on the ModelDB website, with 
accession number 144993. In all the numerical simulations of 
the full stochastic Biophysical neuron model we used Equations 
(1-3) in main text. We used first order Euler-Maruyama inte- 
gration with a time step of dt = 5 /is (quantitative results were 
verified also at dt = 0.5 /xs). Each stimulation pulse was given as 
a square pulse with a width of f s tim = 0.5 ms and amplitude 7o 
(which were respectively named fn and Iq in Soudry and Meir, 





/ [Hz] 

FIGURE 8 | Estimation noise in the cross-power spectral density. To 

estimate the level of this noise in Figure 3D, we added S y j (f) where 
j 7" m J is a shuffled version of { T m \. Only when the estimated Syr (f) is 
above S Y j (f), is its estimation valid. Therefore, in Figure 3D we show only 
this region (left of dashed black line), where estimation is valid. 
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and = 50 ms. Notice the form of the models is given in the 
(more popular) compressed formalism (section 4.1.1), which 
employs the normalization of state occupation probability to 
reduce the dimensionality of equations of Equations (2, 3) in the 
main text. 

4.5.1. The HHS model 

The HHS model combines the Hodgkin-Huxley equations 
(Hodgkin and Huxley, 1952) with slow sodium inactivation 
(Chandler and Meves, 1970; Fleidervish et al, 1996). The model 
equations (Soudry and Meir, 2012b), which employ the uncou- 
pled stochastic noise approximation, are 

CV = g Na sm 3 h (E Na -V)+ g K n 4 (E K - V) 
+gL(E L -V)+I(t) 
m = 4> [a m (V) (1 - m) - fi m (V) m] 

+^N m V (or m (V) (1 - m) + p m (V) m)f m 
h = 4>[a n (V) (1 - n) - p„ (V) n] 

+^N m l 4> («„ (V) (1 - n) + J h (V) n)Hn 
h = 4>[u h (V) (1-h)- /3 h (V) h] 

+^N- 1 4> (a h (V) (l-h) + J h (V) K£h 

s = 5(V)(1 - 5) - y (V) s + ^JVr 1 (8 (V) (l—s) + y (V) s)fe. 

Most of the parameters are given their original values (as in 
Hodgkin and Huxley, 1952; Fleidervish et al., 1996): 

V No = 50mV, V K = -77mV, V L = -54 mV, 
g Na = 120 {kQ ■ cm 2 )' 1 , g K = 36(kQ ■ cm 2 )- 1 , g L = 0.3 {kQ, ■ cm 2 )' 1 , 

c n (V) = 1 °-°^+ 5 ? 55) kHz, MV) = 0.125 ■ e -™/«o kHz, 
«»(V) = £W$m kHz, p m (V) = 4 • e "™/i8 kHZ; 
a h (V) = 0.07 • e -( v+65 V 2() kHz, p h (V) = ( e -°-HV+35) + l)" 1 kHz, 

where in all the rate functions V is used in units of mV. In order to 
obtain the specific spike shape and the latency transients observed 
in cortical neurons, some of the parameters were modified to 

C m = 0.5 /iF/cm 2 , 0 = 2 
Y (V) = 0.51 .L-OMV+17) + A" 1 Hz _ S(V) = 0 .05e-( v+85 >/ 30 Hz. (125) 

We emphasize that these specific choices do not affect any of 
our general arguments, but were chosen for consistency with 
experimental results (Gal et al., 2010). Estimates of channel num- 
ber vary greatly (Soudry and Meir, 2012b). For simplicity, we 
chose N = N n = Nh = N m = N s , and unless stated otherwise, we 
chose, by default N = 10 6 , as in Soudry and Meir (2012b). Note 
that the HHS model is the same model presented in the paper 
with M = 1, </> s ,i = 1> N s ,i = N, N r j = N and <p r = <p. 



4.5.2. The coupled HHS model 

The coupled version of the HHS model uses the same parameters 
as the uncoupled version, and a similar voltage equation 

CV = gNaSomoho (E Na - V) + g K n 0 {E K - V) 
+g L (E L -V)+I(t) 

where the variables «o and sotniho describe the respective frac- 
tion of potassium and sodium channels residing in the "open" 
state. To obtain the coupled model equations, we need to assume 
something about the structure of the ion channels. The original 
assumption by Hodgkin and Huxley was that the channel sub- 
units (e.g., m, n and h) are independent. Over the years, it became 
apparent that this assumption is inaccurate, and the sodium chan- 
nel kinetic subunits are, in fact, not independent (Ulbricht, 2005). 
However, it is not yet clear how the slow sodium inactivation is 
coupled to the rapid channel kinetics (e.g., Menon et al., 2009; 
Milescu et al., 2010), so we nevertheless used the original naive 
HH model assumption that the subunits are independent. In that 
case the potassium channel structure is given by (for brevity, the 
voltage dependence on the rates is henceforth ignored for this 
model) 

4a n 3a n 2a n a n 
n$ =± n\ =± ti2 ^ nj, =± M4 
Pn 2j8„ 3fi n 4/3„ 

while for the sodium channel it is described by 





3 o' jjj ^ c£ 




s 0 m 0 h 0 


^± somiho ^± Sotn 2 ho 


^ som 3 h 0 






SPm 


UhUPh 






a m 


som 0 hi 


=± sotnihi ^± s 0 m 2 hi 


=± s 0 m 3 hi 




Pm ^-Pm 


3Pm 




SllY 






3o; m 2of m 


Urn 


sim 0 h 0 


^ s\tn\ho ^ sim 2 h 0 


^ S\m-$ho 




Pm ^-Pm 


3Pm 


ah 11 Ph 






01 fn 


sitn 0 hi 


=± sirriihi =± sitmhi 


=± sim 3 hi 




Pm ^-Pm 


3Pm 



In this diagram, transition rates indicated between two boxed 
regions, imply that the same rates are used between all corre- 
sponding states in boxed regions. The corresponding 32 SDEs are 
derived using the method described in Orio and Soudry (2012) 
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(or 30 equations if we use the compressed formalism). In this 
model we used Iq = 8.3 /xA. 

4.5.3. The HHSTM model 

In order to investigate the effect of a more "physiological" stimu- 
lation, we changed the HHS model and added synapses. We used 
the popular Tsodyks-Markram model for the effect of a synapse 
with short-term-depression on the somatic voltage (the model 
first appeared in Tsodyks and Markram (1997) and was slightly 
corrected in Tsodyks et al. (1998)). In the model x, y and z are the 
fractions of resources in the recovered, active and inactive states 
respectively, interacting through the system 

y 

/ \ ■ (126) 

x < — z 

Here the z —> x rate is r" 1 , the x — > y rate is r-T , and the x —> 
y rate is UseS (t — t sp ), where S (•) is the Dirac delta function, 
and f sp is the pre-synaptic spike arrival time. The post-synaptic 
current is given by I s (t) = AsEy (t) where Ase is a parameter. 
Additionally, we added noise to the model using the coupled SDE 
method (Orio and Soudry, 2012), assuming that the diagram in 
Equation (126), with the corresponding rates, hint at the under- 
lying Markov kinetic structure, with N = 10 6 . As in Figure IB of 
Tsodyks and Markram (1997), we used x- ln = 3 ms,r rec — 800 ms 
and Use = 0.67. Additionally, we set Ase = 70 fiA to obtain an 
AP response in our model. 

4.5.4. The HHMS model 

The HHMS model consists of many sodium currents, each with 
a different slow kinetic variable. The equations are identical to 
the HHS model, except that g^ a s is replaced by gNaM^ 1 s i> 
where si has the same equation as s in the HHS model, and for 
i > 2, 

5,= [i(V)(l-Si)-y(V)si]^ 

+J(8 (V) (1 - Si ) + y (V) s,) S J <!>,,£, i , 

with 0s, i = £ ' an d N s j = Noe" 1 , where y ( V) and S (V) are taken 
from the HHS model. Unless mentioned otherwise, we chose as 
default 6 = 0.2, r\ = 1.5, M = 5, and N 0 = N as in Figure 5. 

4.5.5. The multiplicative HHMS model 

The Multiplicative HHMS model is identical to the HHMS 
model with r\ = 1, except that gNaM~ l YliL l s ' ' s re pl ace d with 

gNa[[, = iS,. 

4.5.6. The HHSIP model 

The HHSIP model equations (Soudry and Meir, 2012b) are iden- 
tical to the HHS model equations, except that s is renamed to s\ 
and an Inactivating Potassium current was added to the voltage 
equation, where 

JK = gM« 4 s 2 (E K - V) , 
withgvf = 0.0 5gK and 



s 2 = S 2 (V)(l-s 2 )-y 2 {V)s 2 

+Jn~ 1 (8 (V) (1 -s 2 ) + y (V) s 2 )£ s , 2 , 
where N 52 = N and 

3 3e (l/+35)/15 , g -(V+35)/20 

^ W = 1 + <r (v + 35)/io Hz < Yi (V) 

3 3f ,(V+35)/15 + g -(V+35)/20 
1 _|_ e (V+35)/10 

Again, in all the rate functions V is used in mV units. In this 
model we used Iq = 8.3 fiA and T* = 33 ms. 

4.5.7. The HHMSIP model 

The HHMSIP model combines HHSIP and HHMS. Its equations 
are identical to the HHMS model with r) = 2, except they also 
contain the Ik current from the HHSIP model. In this model we 
used Iq = 8.3 fiA and T* = 33 ms, unless otherwise specified. 

ACKNOWLEDGMENTS 

The authors are grateful to O. Barak, N. Brenner, Y. Elhanati, A. 
Gal, T. Knafo, Y. Kafri, S. Marom, and J. Schiller for insightful dis- 
cussions and for reviewing parts of this manuscript. The research 
was partially funded by the Technion V.P.R. fund and by the Intel 
Collaborative Research Institute for Computational Intelligence 
(ICRI-CI). 

REFERENCES 

Anderson, B. D. O., and Moore, J. B. (1979). Optimal Filtering, Vol. 11. Englewood 

Cliffs, NJ: Prentice Hail. 
Bean, B. P. (2007). The action potential in mammalian central neurons. Nat. Rev. 

Neurosci. 8, 451-465. doi: 10.1038/nrn2148 
Bendat, J. S. and Piersol, A. G. (2000). Random Data Analysis and Measurement 

Procedures, Vol. 11, 3rd edn. New York, NY: Wiley. 
Beran, J. (1992). A goodness -of- fit test for time series with long range dependence. 

/. R. Stat. Soc. Ser. B 54, 749-760. 
Beran, J. (1994). Statistics for Long-Memory Processes. New York, NY: Chapman & 

Hall. 

Brecht, M., Schneider, M., Sakmann, B., and Margrie, T. W. (2004). Whisker move- 
ments evoked by stimulation of single pyramidal cells in rat motor cortex. 
Nature 427, 704-710. doi: 10.1038/nature02266 

Chandler, W. K., and Meves, H. (1970). Slow changes in membrane permeability 
and long-lasting action potentials in axons perfused with fluoride solutions. /. 
Physiol 211,707-728. 

Channelpedia. Available online at: http://channelpedia.epfl.ch/ 

Colquhoun, D., and Hawkes, A. G. (1981). On the stochastic properties of sin- 
gle ion channels. Proc. R. Soc. Lond. Ser. B Biol Sci. 211, 205-235. doi: 
10.1098/rspb.l981.0003 

Contou-Carrere, M. N. (2011). Model reduction of multi-scale chemical Langevin 
equations. Syst. Control Lett. 60, 75-86. doi: 10.1016/j.sysconle.2010.10.011 

De Col, R., Messlinger, K., and Carr, R. W. (2008). Conduction velocity is regu- 
lated by sodium channel inactivation in unmyelinated axons innervating the rat 
cranial meninges. 7. Physiol 586, 1089-1103. doi: 10.1113/jphysiol.2007.145383 

De Paola, V., Holtmaat, A, Knott, G., Song, S., Wilbrecht, L., Caroni, P., et al. 
(2006). Cell type-specific structural plasticity of axonal branches and boutons 
in the adult neocortex. Neuron 49, 861-875. doi: 10.1016/j.neuron.2006.02.017 

Debanne, D., Campanac, E., Bialowas, A, and Carlier, E. (201 1). Axon physiology. 
Physiol Rev. 91, 555-602. doi: 10.1152/physrev.00048.2009 

Druckmann, S., Berger, T. K., Schurmann, E, Hill, S., Markram, H., and Segev, I. 
(201 1). Effective stimuli for constructing reliable neuron models. PLoS Comput. 
Biol. 7:el002133. doi: 10.1371/journal.pcbi.l002133 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



April 2014 | Volume 8 | Article 29 | 21 



Soudry and Meir 



A linearized spiking input-output 



Elul, R., and Adey, W. R. (1966). Instability of firing threshold and "remote" 

activation in cortical neurons. Nature 212, 1424-1425. doi: 10.1038/2121424a0 
Ermentrout, B., and Terman, D. (2010). Mathematical Foundations of Neuroscience, 

Vol. 35. New York, NY: Springer, doi: 10.1007/978-0-387-87708-2 
Fleidervish, I. A., Friedman, A., and Gutnick, M. J. (1996). Slow inactivation of 

Na+ current and slow cumulative spike adaptation in mouse and guinea-pig 

neocortical neurones in slices. /. Physiol. 493, 83-97. 
Fox, R. F., and Lu, Y. N. (1994). Emergent collective behavior in large numbers 

of globally coupled independently stochastic ion channels. Phys. Rev. E 49, 

3421-3431. doi: 10.1103/PhysRevE.49.3421 
Gal, A., Eytan, D., Wallach, A., Sandler, M., Schiller, J., and Marom, S. (2010). 

Dynamics of excitability over extended timescales in cultured cortical neurons. 

/. Neurosci. 30, 16332-16342. doi: 10.1523/JNEUROSCI.4859-10.2010 
Gal, A., and Marom, S. (2013). Entrainment of the intrinsic dynamics of sin- 
gle isolated neurons by natural-like input. /. Neurosci. 33, 7912-7918. doi: 

10.1523/JNEUROSCI.3763-12.2013 
Gardiner, C. W. (2004). Handbook of Stochastic Methods, 3rd edn. Berlin: Springer- 

Verlag. doi: 10.1007/978-3-662-05389-8 
Gerstner, W., and Naud, R. (2009). How good are neuron models? Science 326, 

379-380. doi: 10.1126/science.ll81936 
Goldwyn, J. H., Imennov, N. S., Famulare, M., and Shea-Brown, E. (2011). 

Stochastic differential equation models for ion channel noise in Hodgkin- 

Huxley neurons. Phys. Rev. E 83, 041908. doi: 10.1 103/PhysRevE.83.041908 
Goldwyn, J. H., Rubinstein, J. X, and Shea-Brown, E. (2012). A 

point process framework for modeling electrical stimulation of the 

auditory nerve. /. Neurophysiol. 108, 1430-1452. doi: 10.1152/jn. 

00095.2012 

Grubb, M. S., and Burrone, J. (2010). Activity-dependent relocation of the axon 
initial segment fine-tunes neuronal excitability. Nature 465, 1070-1074. doi: 
10.1038/nature09160 

Hille, B. (2001). Ion Channels of Excitable Membranes, 3rd edn. Sunderland, MA: 
Sinauer Associates. 

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane 
current and its application to conduction and excitation in nerve. /. Physiol. 
117, 500. 

Huys, Q. J. M., Ahrens, M. B., and Paninski, L. (2006). Efficient esti- 
mation of detailed single-neuron models. /. Neurophysiol. 96, 872. doi: 
10.1152/jn.00079.2006 

Ikegaya, Y, Sasaki, T., Ishikawa, D., Honma, N., Tao, K., Takahashi, N., 
et al. (2012). Interpyramid spike transmission stabilizes the sparseness of 
recurrent network activity. Cereb. Cortex 23, 293-304. doi: 10.1093/cercor/ 
bhs006 

Jugloff, D. G. M. (2000). Internalization of the Kvl.4 potassium channel is sup- 
pressed by clustering interactions with PSD-95. /. Biol. Chem. 275, 1357-1364. 
doi: 10.1074/jbc.275.2.1357 

Kaplan, D. T., Clay, J. R., Manning, T., Glass, L., Guevara, M. R., and 
Shrier, A. (1996). Subthreshold dynamics in periodically stimulated squid 
giant axons. Phys. Rev. Lett. 76, 4074-4077. doi: 10.1103/PhysRevLett. 
76.4074 

Kasischke, K. A., Vishwasrao, H. D. D., Fisher, P. J., Zipfel, W. R., and Webb, 
W. W. (2004). Neural activity triggers neuronal oxidative metabolism followed 
by astrocytic glycolysis. Science (New York, NY.) 305, 99-103. doi: 10.1126/sci- 
ence. 1096485 

Keshner, M. S. (1982). 1/f noise. Proc. IEEE 70, 212-218. doi: 
10.1109/PROC.1982.12282 

Koch, C, and Segev, I. (1989). Methods in Neuronal Modeling: From Ions to 
Networks, Vol. 484, 2nd edn. Cambridge: MIT Press. 

Komlosi, G., Molnar, G., Rozsa, M., Olah, S., Barzo, P., and Tamas, G. 
(2012). Fluoxetine (prozac) and serotonin act on excitatory synaptic trans- 
mission to suppress single layer 2/3 pyramidal neuron-triggered cell assem- 
blies in the human prefrontal cortex. /. Neurosci. 32, 16369-16378. doi: 
10.1523/JNEUROSCI.2618-12.2012 

Lejung, L. (1999). System Identification: Theory for the User. Upper Saddle River, 
NJ: PTR Prentice Hall. 

Levitan, I. B. (1994). Modulation of ion channels by protein phospho- 
rylation and dephosphorylation. Annu. Rev. Physiol. 11, 193-212. doi: 
10.1146/annurev.ph.56.030194.001205 

Li, C. Y. T., Poo, M. M., and Y D. (2009). Burst spiking of a single cortical neuron 
modifies global brain state. Science 324, 643-646. doi: 10. 1126/science.l 169957 



Linaro, D., Storace, M., and Mattia, M. (2011). Inferring network dynamics and 

neuron properties from population recordings. Front. Comput. Neurosci. 5:43. 

doi: 10.3389/fncom.201 1.00043 
Marom, S. (2010). Neural timescales or lack thereof. Prog. Neurobiol 90, 16-28. 

doi: 10.1016/j.pneurobio.2009.10.003 
Menon, V., Spruston, N., and Kath, W. L. (2009). A state-mutating genetic 

algorithm to design ion-channel models. Proc. Natl. Acad. Sci. U.S.A. 106, 

16829-16834. doi: 10.1073/pnas.0903766106 
Migliore, M., Cannia, C, Lytton, W. W., Markram, H., and Hines, M. L. (2006). 

Parallel network simulations with NEURON. /. Comput. Neurosci. 21,1 19-129. 

doi: 10.1007/sl0827-006-7949-5 
Milescu, L. S., Yamanishi, T., Ptak, K., and Smith, J. C. (2010). Kinetic 

properties and functional dynamics of sodium channels during repetitive 

spiking in a slow pacemaker neuron. /. Neurosci. 30, 12113-12127. doi: 

10.1523/JNEUROSCI.0445-10.2010 
Molnar, G., Olah, S., Komlosi, G., Fiile, M., Szabadics, J., Varga, C, et al. (2008). 

Complex events initiated by individual spikes in the human cerebral cortex. 

PLoSBiol. 6:e222. doi: 10.1371/journal.pbio.0060222 
Monjaraz, E., Navarrete, A., Lopez-Santiago, L. E, Vega, A. V, and Cota, G. (2000). 

L-type calcium channel activity regulates sodium channel levels in rat pituitary 

GH3 cells. /. Physiol. 523, 45-55. doi: 10.1111/j.l469-7793.2000.00045.x 
Neher, E., and Sakmann, B. (1976). Single-channel currents recorded from 

membrane of denervated frog muscle fibres. Nature 260, 799-802. doi: 

10.1038/260799a0 

Nishiyama, H., Fukaya, M., and Watanabe, M. (2007). Axonal motil- 
ity and its modulation by activity are branch-type specific in the 
intact adult cerebellum. Neuron 56, 472-487. doi: 10. 1016/j. neuron. 
2007.09.010 

Orio, P., and Soudry, D. (2012). Simple, fast and accurate implementation of the 
diffusion approximation algorithm for stochastic ion channels with multiple 
states. PLoS ONE 7:e36670. doi: 10.1371/journal.pone.0036670 

Papoulis, A., and Pillai, S. U. (1965). Probability, Random Variables, and Stochastic 
Processes. New York, NY: McGraw-Hill. 

Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applica- 
tions in speech recognition. Proc. IEEE 77, 257-286. doi: 10.1 109/5.18626 

Roth, A., and Hausser, M. (2001). Compartmental models of rat cerebel- 
lar Purkinje cells based on simultaneous somatic and dendritic patch- 
clamp recordings. /. Physiol. 535, 445-472. doi: 10.111 1469-7793.2001. 
00445.x 

Schneidman, E., Freedman, B., and Segev, I. (1998). Ion channel stochas- 
ticity may be critical in determining the reliability and precision of 
spike timing. Neural Comput. 10, 1679-1703. doi: 10.1162/089976698 
300017089 

Silver, I. A., Deas, J., and Erecinska, M. (1997). Ion homeostasis in brain cells: 
differences in intracellular ion responses to energy limitation between cul- 
tured neurons and glial cells. Neuroscience 78, 589-601. doi: 10.1016/S0306- 
4522(96)00600-8 

Sjostrom, P. J. J., Rancz, E. A. A., Roth, A., and Hausser, M. (2008). Dendritic 
excitability and synaptic plasticity. Physiol. Rev. 88, 769-840. doi: 10.1 152/phys- 
rev.00016.2007 

Song, S., Sjostrom, P. J. J., Reigl, M., Nelson, S. B., and Chklovskii, D. B. (2005). 

Highly nonrandom features of synaptic connectivity in local cortical circuits. 

PLoSBiol 3:e68. doi: 10.1371/journal.pbio.0030068 
Soudry, D., and Meir, R. (2012a). An exact reduction of the master equation to a 

strictly stable system with an explicit expression for the stationary distribution. 

arXiv: 1207.4436. 

Soudry, D., and Meir, R. (2012b). Conductance-based neuron models and 
the slow dynamics of excitability. Front. Comput. Neurosci. 6:4. doi: 
10.3389/fncom.2012.00004 

Soudry, D., and Meir, R. (2014). The neuronal response at extended timescales: 
long term correlations without long memory. Front. Comput. Neurosci. 8:35. 
doi: 10.3389/fncom.2014.00035 

Staub, O., Gautschi, I., Ishikawa, X, Breitschopf, K., Ciechanover, A., Schild, 
L., et al. (1997). Regulation of stability and function of the epithelial 
Na+ channel (ENaC) by ubiquitination. EMBO J. 16, 6325-6336. doi: 
10.1093/emboj/16.21.6325 

Toib, A., Lyakhov, V., and Marom, S. (1998). Interaction between duration of activ- 
ity and time course of recovery from slow inactivation in mammalian brain Na+ 
channels. /. Neurosci. 18, 1893-1903. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



April 2014 | Volume 8 | Article 29 | 22 



Soudry and Meir 



A linearized spiking input-output 



Tsodyks, M., and Markram, H. (1997). The neural code between neo- 
cortical pyramidal neurons depends on neurotransmitter release prob- 
ability. Proc. Natl. Acad. Sri. U.S.A. 94, 719-723. doi: 10.1073/pnas. 
94.2.719 

Tsodyks, M., Pawelzik, K., and Markram, H. (1998). Neural networks with 
dynamic synapses. Neural Comput. 10, 821-835. doi: 10.1162/08997669830 
0017502 

Ulbricht, W. (2005). Sodium channel inactivation: molecular determinants and 
modulation. Physiol Rev. 85, 1271-1301. doi: 10.1152/physrev.00024.2004 

Wainrib, G., Thieullen, M., and Pakdaman, K. (2011). Reduction of 
stochastic conductance-based neuron models with time-scales sep- 
aration. /. Comput. Neurosci. 32, 327-346. doi: 10.1007/sl0827- 
011-0355-7 

Wallach, A. (2012). The Response Clamp : A Control Based Approach for 
the Study of Neural Systems; Method and Applications. Ph.D. thesis, 
Technion. 



Conflict of Interest Statement: The authors declare that the research was con- 
ducted in the absence of any commercial or financial relationships that could be 
construed as a potential conflict of interest. 

Received: 20 December 2013; accepted: 24 February 2014; published online: 02 April 
2014. 

Citation: Soudry D and Meir R (2014) The neuronal response at extended timescales: a 
linearized spiking input-output relation. Front. Comput. Neurosci. 8:29. doi: 10.3389/ 
fncom.2014.00029 

This article was submitted to the journal Frontiers in Computational Neuroscience. 
Copyright © 2014 Soudry and Meir. This is an open-access article distributed under 
the terms of the Creative Commons Attribution License (CC BY). The use, distribu- 
tion or reproduction in other forums is permitted, provided the original author (s) 
or licensor are credited and that the original publication in this journal is cited, in 
accordance with accepted academic practice. No use, distribution or reproduction is 
permitted which does not comply with these terms. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



April 2014 | Volume 8 | Article 29 | 23 



